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Abstract 

A detailed treatment of the classical Chapman-Enskog derivation 
of hydrodynamics is given in the framework of Grad's moment equa- 
tions. Grad's systems are considered as the minimal kinetic models 
where the Chapman-Enskog method can be studied exactly, thereby 
providing the basis to compare various approximations in extending 
the hydrodynamic description beyond the Navier-Stokes approxima- 
tion. Various techniques, such as the method of partial summation, 
Pade approximants, and invariance principle are compared both in lin- 
ear and nonlinear situations. 



* Corresponding author. Email: ikarlin@ifp.mat.ethz.ch Fax: +41 1 632 10 76 
^Work supported by the Clay Mathematics Institute 



1 



Contents 

1 Introduction 

1.1 The 'ultra-violet catastrophe' of the Chapman-Enskog expansion 

1.2 The Chapman-Enskog method for linearized Grad's equations 

2 Exact summation of the Chapman-Enskog expansion 

2.1 The IDIOM Grad's equations 

2.2 The 3D10M Grad's equations 

3 The dynamic invariance principle 

3.1 Partial summation of the Chapman-Enskog expansion 

3.2 The dynamic invariance 

3.3 The Newton method 

3.4 The invariance equation for the 1D13M Grad system 

3.5 The invariance equation for the 3D13M Grad system 

3.6 Gradient expansion in the kinetic theory of phonons 

3.7 Nonlinear Grad equations. 

4 Conclusion 
References 



2 



1 Introduction 



1.1 The "ultra-violet catastrophe" of the Chapman-Enskog 
expansion 

Most of the interesting expansions in non-equilibrium statistical physics are 
divergent. This paraphrase of the well known folklore "Dorfman's theorem" 
conveys the intrinsic problem of many-body systems: A number of system- 
atic (at the first glance) methods has led to (a). An excellent but already 
known on the phenomenological grounds first approximation, (b). Already 
the next correction, not known phenomenologically and hence of interest, 
does not exist because of a divergence. There are many examples of this 
situations: Cluster expansion of the exact collision integral for dense gases 
lead to divergent approximations of transport coefficients, non-convergent 
long tails of correlation functions in the Green-Kubo formulae etc. 

Derivation of the hydrodynamic equations from a microscopic description 
is the classical problem of physical kinetics. As is well known, the famous 
Chapman-Enskog method [1] provides an opportunity to compute a solution 
from the Boltzmann kinetic equation as a formal series in powers of the 
Knudsen number e. The parameter e reflects the ratio between the mean 
free path of a particle, and the scale of variations of the hydrodynamic 
fields, the density, the mean flux, and the temperature. If the Chapman- 
Enskog expansion is truncated at a certain order, we obtain subsequently: 
the Euler hydrodynamics ( e°), the Navier-Stokes hydrodynamics (e 1 ), the 
Burnett hydrodynamics (e 2 ), the super-Burnett hydrodynamics (e 3 ), etc. 
The post-Navier-Stokes terms extend the hydrodynamic description beyond 
the strictly hydrodynamic limit eCl. 

However, as it has been first demonstrated by Bobylev [2], even in the 
simplest regime (one-dimensional linear deviations around the global equilib- 
rium), the Burnett hydrodynamic equations violate the basic physics behind 
the Boltzmann equation. Namely, sufficiently short acoustic waves are in- 
creasing with time instead of decaying. This contradicts the .ff-theorem, 
since all near-equilibrium perturbations must decay. The situation does not 
improve in the next, super-Burnett approximation. 

This "ultra-violet catastrophe" which occurs in the lower-order trunca- 
tions of the Chapman-Enskog expansion creates therefore very serious diffi- 
culties in the problem of an extension of the hydrodynamic description into 
a highly non-equilibrium domain (see [3] for a discussion of other difficulties 
of the post-Navier-Stokes terms of the Chapman-Enskog expansion). The 
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Euler and the Navier-Stokes approximations remain basic in the hydrody- 
namic description, while the problem of their extension is one of the central 
open problems of the kinetic theory. The study of approximate solutions 
based on the Chapman-Enskog method still continues [4]. 

All this begs for a question: What is wrong with the Chapman-Enskog 
method? At the first glance, the failure of the Burnett and of the super- 
Burnett hydrodynamics may be accounted in favor of a frequently used 
argumentation about the asymptotic character of the Chapman-Enskog ex- 
pansion. However, it is worthwhile to notice here that divergences in the 
low-order terms of formal expansions are not too surprising. In many oc- 
casions, in particular, in quantum field theory [5] and in statistical physics 
[6] , the situation is often improved if one takes into account the very remote 
terms of corresponding expansions. Thus, a more constructive viewpoint 
on the Chapman-Enskog expansion could be to proceed along these lines, 
and to try to sum up the Chapman-Enskog series, at least formally and 
approximately. 

An attempt of this kind of working with the Chapman-Enskog expansion 
is undertaken in this paper. The formalities are known to be rather awk- 
ward for the Boltzmann equation, and untill now, exact summations of the 
Chapman-Enskog expansion are known in a very limited number of cases 
[7]. In this paper, we will concentrate on the Chapman-Enskog method as 
applied to the well known Grad moment equations [8] . The use of the Grad 
equations for our purpose brings, of course, considerable technical simplifi- 
cations as compared to the case of the Boltzmann equation but it does not 
make the problem trivial. Indeed, the Chapman-Enskog method amounts to 
a nonlinear recurrence procedure even as applied to the simplest, linearized 
Grad equations. Moreover, as we will see it soon, the Chapman-Enskog ex- 
pansion for moment systems inherits Bobylev's instability in the low-order 
approximations. Still, the advantage of our approach is that many explicit 
results can be obtained and analyzed. In order to summarize, in this pa- 
per we consider Grad's moment equations as finitely-coupled kinetic mod- 
els where the problem of reduced description is meaningful, rather than as 
models of extended hydrodynamics. The latter viewpoint is well known as a 
microscopic background of the extended irreversible thermodynamics [9, 10]. 

The outline of this paper is as follows: after an introduction of the 
Chapman-Enskog procedure for the linearized Grad equations (section 1.2), 
we will open the discussion with two examples (the linearized one- and three- 
dimensional 10 moment Grad equations) where the Chapman-Enskog series 
is summed up exactly in a closed form (sections 2.1 and 2.2). These results 
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makes it possible to discuss the features of the Chapman-Enskog solution 
in the short-wave domain in the frames of the model, and will serve for a 
purpose of testing various approximate methods thereafter. We will see, in 
particular, that the "smallness" of the Knudsen number e used to develop 
the Chapman-Enskog method has no direct meaning in the exact result. 
Also, it will become clear that the finite-order truncations, even provided 
they are stable, give less opportunities to approximate the solution in a 
whole, and especially in the short-wave domain. 

The exact solutions are, of course, the lucky exceptions, and even for the 
Grad moment equations the complexity of the Chapman-Enskog method 
increases rapidly with an increase of the number of the moments taken into 
account. Further (section 3.1) we will review a technique of summing the 
Chapman-Enskog expansion partially. This technique is heuristic (as are 
the methods of partial summing in general) , but still it removes the Bobylev 
instability, as well as it qualitatively reproduces the features of the exact 
solutions in the short-wave limit. 

The style of working in the sections mentioned so far falls into a paradigm 
of the Taylor- like expansions into powers of the Knudsen number. This view- 
point on the problem of a derivation of the hydrodynamics will be altered 
beginning with the section 3.2. There we demonstrate that a condition of a 
dynamic invariance which can be realized directly and with no request of the 
Knudsen number brings us to the same result as the exact summation of the 
Chapman-Enskog expansion. The Chapman-Enskog method thereafter can 
be regarded for a one possibility to solve the resulting invariance equations. 
Further, we demonstrate that iterative methods provide a reasonable alter- 
native to the Taylor expansion in this problem. Namely, we show that the 
Newton method has certain advantages above the Chapman-Enskog method 
(section 3.3). We also establish a relationship between the method of partial 
summing and the Newton method. 

The material of further sections serves for an illustrative introduction 
how the pair 'invariance equation + Newton method' can be applied to the 
problems of kinetic theory. The remaining sections of this chapter are de- 
voted to further examples of this type of working on the level of the Grad 
equations. In sections 3.4 and 3.5 we derive and discuss the invariance 
equations for the linearized thirteen-moment Grad equations. Section 3.6 
is devoted to kinetic equations of the Grad type, arising in the problems 
of phonon transport in massive solids at low temperatures. In particular, 
we demonstrate that the onset of the second sound regime of phonon prop- 
agation corresponds to a branching point of the exact sum of the relevant 
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Chapman-Enskog expansion. 

In section 3.7 we apply the invariance principle to nonlinear Grad equa- 
tions. We sum up exactly a subseries of the Chapman-Enskog expansion, 
namely, the dominant contribution in the limit of high average velocities. 
This type of contribution is therefore important for an extension of the hy- 
drodynamic description into the domain of strong shock waves. We give a 
relevant analysis of the corresponding invariance equation, and, in particu- 
lar, discuss the nature of singular points of this equation. A brief discussion 
concludes this paper. Some of the results presented below were published 
earlier in [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] 

1.2 The Chapman-Enskog method for linearized Grad's equa- 
tions 

In this section, for the sake of completeness, we introduce linearized Grad's 
equations and the Chapman-Enskog method for them in a form to be used in 
this paper. Since the Chapman-Enskog method is extensively discussed in a 
number of books, especially, in the classical monograph [1], our presentation 
will be brief. 

Notations will follow those of the papers [2, 11]. We denote po, Tq and 
u = as the fixed equilibrium values of the density, of the temperature and 
of the averaged velocity (in the appropriate Galilean reference frame) , while 
5p, 5T and Su are small deviations of the hydrodynamic quantities from 
their equilibrium values. Grad's moment equations [8] which will appear 
below, contain the temperature dependent viscosity coefficient p(T). It is 
convenient to write p(T) = rj(T)T. The functional form of t](T) is dictated 
by the choice of the model for particle's interaction. In particular, we have 
rj = const for Maxwell's molecules, and rj ~ VT for hard spheres. Let us 
introduce the following system of dimensionless variables: 

Su Sp ST 

u = , p = — , T=—, (1) 

^k B T /m po Jo 

x = P° t=-^t' 

V (T )^k B T /m ' r)(T ) ' 

where x' are spatial coordinates, and t' is the time. In the sequel, we use the 
system of units in which Boltzmann's constant k B and the particle's mass 
m are one. Three-dimensional thirteen moment Grad's equations, linearized 
near the equilibrium, take the following form when written in terms of the 
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dimensionless variables (1): 



dtp = 


-V • u, 


(2) 


d t u = 


-Vp- VT- V • <x, 




d t T = 


-^(V-ti + V-g), 




d t <J = 


2 

-Vu - -Vq — rr. 

5 


(3) 


d t q = 


--VT- V-er- -g. 
2 3 





In these equations, cr(x,t) and q(x,t) are dimensionless quantities corre- 
sponding to the stress tensor and to the heat flux, respectively. Further, the 
gradient V stands for the vector of spatial derivatives d/dx. The dot de- 
notes the standard scalar product, while the overline stands for a symmetric 
traceless dyad. In particular, 

2 

Vu = Vu + (Vu) T - -IV ■ u, 

where / is unit matrix. 

Grad's equations (2) and (3) is the simplest model of a coupling of the hy- 
drodynamic variables, p(x,t), T(x,t) and u(x,t), to the non-hydrodynamic 
variables cr(x,t) and q(x,t). The problem of reduced description is to close 
the first three equations (2), and to get an autonomous system for the hydro- 
dynamic variables alone. In other words, the non-hydrodynamic variables 
tr(x,t) and q(x,t) should be expressed in terms of the variables p(x,t), 
T(x,t) and u(x,t). The Chapman-Enskog method, as applied for this pur- 
pose to Grad's system (2) and (3), involves the following steps: 

First, we introduce a formal parameter e, and write instead of equations 

(3): 

2 1 

d t cr = -Vu - -Vq - -a, (4) 

O 6 

t q = — -VT — V • <r — —q. 
H 2 3e 

Second, the Chapman-Enskog solution is found as a formal expansions of 

the stress tensor and of the heat flux vector: 

oo 

a = ^e n+ V n ); (5) 

n=0 

oo 

q = £ e "+y«). 

n=0 
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Zero-order coefficients, cr^ and q^ are: 



<r(°)=-W, g(°) = -^VT. (6) 

Coefficients of the order n > 1 are found from the recurrence procedure: 

r n-l 



(n) 



>) = A 

' 2 



^ Q^^in-l-m) + f V g(n-1) j (?) 
m=0 5 J 

^q^-^+V-o-^- 1 )}, 

lm=0 J 



where c^" 1 ^ are recurrently defined Chapman- Enskog operators. They act 
on functions p(x,t), T(x,t) and u(x,t), and on their spatial derivatives, 
according to the following rule: 

<T> D , :>?; (B) 

M _ f -§Z)V-« m = . 

dt U1 ~ \ -lDV-q^-V m>\ ' 

(m) / -DV(p + T) m = 

d * ^ " \ -DV ■ m > 1 ■ 

Here D is an arbitrary differential operator with constant coefficients. 

Given the initial condition (6), the Chapman-Enskog equations (7) and 
(8) are recurrently solvable. Finally, by terminating the computation at the 
order N > 0, we obtain the iVth order approximations to the expansions 
(5), a N and q N : 

N N 

<^v = E en+1<T(n) ' <^ = E en+ V n) - (9) 

n=0 n=0 

Substituting these expressions instead of the functions er and q in the equa- 
tions (2), we close the latter to give the hydrodynamic equations of the 
order N. In particular, N = results in the Navier-Stokes approximation, 
N = 1 and N = 2 give the Burnett and the super-Burnett approximations, 
respectively, and so on. 

Though the 'microscopic' features of Grad's moment equations are, of 
course, much simpler as compared with the Boltzmann equation, the Chapman- 
Enskog procedure for them just described is not trivial. Our purpose is to 
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study explicitly the features of the gradient expansions like (5) in the highly 
non-equilibrium domain, in particular, to find out to what extend the finite- 
order truncations (9) approximate the solution, and what kind of alternative 
strategies to find approximations are possible. In the following, when refer- 
ring to Grad's equations, we use the notation mDnM, where m is the spatial 
dimension of corresponding fields, and n is the number of these fields. For 
example, the above system is the 3D13M Grad's system. 

2 Exact summation of the Chapman— Enskog ex- 
pansion 

2.1 The IDIOM Grad equations 

In this section, we open our discussion with the exact summation of the 
Chapman-Enskog series for the simplest Grad's system, the one-dimensional 
linearized ten-moment equations. Throughout the section we use the hydro- 
dynamic variables p(x,t) = p(x,t) + T(x,t) and u(x,t), representing the 
dimensionless deviations of the pressure and of the average velocity from 
their equilibrium values. The starting point is the linearized Grad's equa- 
tions for p, u, and a, where a is the dimensionless xx-component of the 
stress tensor: 



The system of equations for three functions is the derivative of the ten- 
moment Grad's system (see Eq. (38) below). Equations (10) provides the 
simplest model of a coupling of the hydrodynamic variables, u and p, to 
the single non-hydrodynamic variable a, and corresponds to a heat non- 
conductive case. 

Our goal here is to shorten the description, and to get a closed set of 
equations with respect to variables p and u only. That is, we have to express 
the function a in the terms of spatial derivatives of the functions p and of u. 
The Chapman-Enskog method, as applied to eq. (10) results in the following 



d t p 
d t u 



d x p ~ d x a, 



(10) 
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series representation: 

oo 

a = ^e n+ V n) . (11) 

n=0 

Coefficients are due to the following recurrence procedure [11]: 

n— 1 

(7 W=-^fi{ m ) (7 (n-l-m) i (12 ) 
m=0 

where the Chapman-Enskog operators act on p, on u, and on their 

spatial derivatives as follows: 

flHfli,. _ / m = r . 

^ " j m>l ' (13) 

1 21 10, m > 1 

Here / > is an arbitrary integer, and <9° = 1. Finally, 

= -^u, (14) 

which leads to the Navier-Stokes approximation of the stress tensor: ctns = 
ea(°). 

Because of a somewhat involved structure of the recurrence procedure 
(12) and (13), the Chapman-Enskog method is a nonlinear operation even in 
the simplest model (10). Moreover, the Bobylev instability is again present. 

Indeed, computing the coefficients crW and on the basis of the ex- 
pressions (12), we obtain: 

a B = e< 7<°> + eV 1 ) = ~{ed x u + e 2 d 2 xP ), (15) 

and 

a SB = e<7 (°) + eV 1 ) + eV 2 ) = -|(e0 x u + e 2 d 2 p + \e 3 d 3 x u), (16) 

for the Burnett and the super-Burnett approximations, respectively. Now 
we can substitute each of the approximations, ctns> °"b, and ctsb instead 
of the function a in the second of the equations of the set (10). Then 
the equations thus obtained, together with the equation for the density 
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p, will give us the closed systems of the hydrodynamic equations of the 
Navier-Stokes, of the Burnett, and of the super-Burnett levels. To see the 
properties of the resulting equations, we compute the dispersion relation for 
the hydrodynamic modes. Using a new space-time scale, x' = e _1 x, and 
t' = e~ 1 t, and next representing u = Uk^p(x',t'), and p = Pk<f(x' ,t'), where 
ip(x',t') = exp((jt' + ikx'), and k is a real-valued wave vector, we obtain 
the following dispersion relations u(k) from the condition of a non-trivial 
solvability of the corresponding linear system with respect to and p^- 

co ± = -h 2 ±h\k\V^-15, (17) 
for the Navier-Stokes approximation, 

u± = --k 2 ±-i\k\V8k 2 + 15, (18) 
for the Burnett approximation (15), and 

u± = ^k 2 (k 2 - 3) ± h\k\ V4:k® - 24k 4 - 72k 2 - 135, (19) 

for the super-Burnett approximation (16). 

These examples demonstrate that the real part Re u±{k) < for the 
Navier-Stokes (17) and for the Burnett (18) approximations, for all wave 
vectors. Thus, these approximations describe attenuating acoustic waves. 
However, for the super-Burnett approximation, the function Re uj±(k) (19) 
becomes positive as soon as \k\ > y/3. That is, the equilibrium point is stable 
within the Navier-Stokes and the Burnett approximation, and it becomes un- 
stable within the super-Burnett approximation for sufficiently short waves. 
Similar to the case of the Bobylev instability of the Burnett hydrodynam- 
ics for the Boltzmann equation, the latter result contradicts the dissipative 
properties of the Grad system (10): the spectrum of the full IDIOM system 
(10) is stable for arbitrary k. 

Our goal now is to sum up the series (11) in a closed form. Firstly, we 
will make some preparations. 

As was demonstrated in [11], the functions in eqs. (11) and (12) 
have the following explicit structure to arbitrary order n > 0: 

= a n d 2n+1 u, (20) 
a( 2 ™ +1 ) = b n d 2 x ^p, 
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where the coefficients a n and b n are determined through the recurrence pro- 
cedure (12), and (13). The Chapman-Enskog procedure (12) and (13) can 
be represented in terms of the real- valued coefficients a n and b n (20). We 
will do this below. 

Knowing the structure (20) of the coefficients of the Chapman-Enskog 
expansion (11), we can write down its formal sum. It is convenient to use the 
Fourier variables introduced above which amounts essentially to the change 
ed x — > ik. Substituting the expression (20) into the Chapman-Enskog series 
(11), we obtain the following formal expression for the Fourier image of the 
sum: 

a k = ikA{k 2 )u k - k 2 B(k 2 ) Pk , (21) 

where the functions A(k 2 ) and B(k 2 ) are formal power series with the coef- 
ficients (20): 

oo 

A(k 2 ) = J2^n(-k 2 ) n , (22) 

n=0 

oo 

B(k 2 ) = Y,U-k 2 ) n . 

n=0 

Thus, the question of the summation of the Chapman-Enskog series (11) 
amounts to finding the two functions, A(k 2 ) and B{k 2 ) (22). Knowing A 
and B, we derive a dispersion relation for the hydrodynamic modes: 



u± = ^±^k*A*-^{l-k*B). (23) 

Now we will concentrate on the problem of a computation of the functions 
A and B (22) in a closed form. For this purpose, we will first express the 
Chapman-Enskog procedure (12) and (13) in terms of the coefficients a n 
and b n (20). At the same time, our derivation will constitute a proof of the 
structure (20). 

It is convenient to start with the Fourier representation of the equa- 
tions (12) and (13). Writing u = u k exp(ifcx), p = pfcexp(ifcx), and a = 
<7fc exp(ifcx), we obtain the following representation: 

9t Uk = j -ik*T X \ rn>l ' (24) 

a (m) f -hku k , m = 

d < Pk = \ 0, m>l ' 
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while 

n-1 

4 n) = -£^ m) 4 n " 1-m) , (25) 

m=0 

and 



4^ = a n (-k 2 ) n iku k , (26) 
+1 ) = 6n (-fe 2 ) n (-fe 2 )p fc . 

The Navier-Stokes and the Burnett approximations give ao = — |, and 
bo = — |. Thus, the structure (26) is proven for n = 0. 

The further derivation relies upon induction. Let us assume that the 
ansatz (26) is proven up to the order n. Computing the coefficient a^ n+1 ^ 
from the equation (25), we have: 



it it 

a (2(n+l)) = _ d (0) (j (2n+l) _ y- Q {2m+1) ^(n-m)) _ y- g ( 
h t h / j t h / j t 



,(2m) (2(n-m)+l) 
m=0 m=l 

(27) 

Due to the assumption of the induction, we can adopt the form of the 
coefficients (26) in all the terms on the right hand side of the latter 
expression. On the basis of the expressions (26) and (24), we conclude that 
each summand in the last sum in (27) is equal to zero. Further, the term 
df^ a^ n+1 ^ gives a linear contribution: 

while the terms in the remaining sum contribute nonlinearly: 

<9 t (2m+1 M 2(n ~ m)) = a n _ m (-A^-^H (2m+1 V = -a n . m a m {-k^Hku k . 

Substituting the last two expressions into the equation (27), we see that it 
has just the same structure as the coefficient a^ 2 ^ 1 ^ in (26). Thus, we 
come to the first recurrence equation: 

5 n 

O-n+l = ~^b n + ^ ' (l n — m (l m . 

6 m=0 

Computing the coefficient a ^ n+1 ^ +1 ^ by the same pattern, we come to the 
second recurrence equation, and the Chapman-Enskog procedure (12) and 
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(13) becomes the following equivalent formulation in terms of the coefficients 
a n and b n (20): 

5 n 

O-n+l = + E a n-md m , (28) 

n 

b n +i = o-n+i + ^ a, n - m b m . 

m=0 

The initial condition for this set of equations is dictated by the Navier-Stokes 
and the Burnett terms: 

4 4 
°o = -3. 6 ° = -3 ( 29 ) 

Our goal now is to compute the functions A and B (22) on the basis of 
the recurrence equations (28). At this point, it is worthwhile to notice that 
usual routes of dealing with the recurrence system (28) would be either to 
truncate it at a certain n, or to calculate all the coefficients explicitly, and 
next to substitute the result into the power series (22). Both these routes 
are not successful here. Indeed, retaining the coefficients ao, bo, and a\ 
gives the super-Burnett approximation (16) which has the Bobylev short- 
wave instability, and there is no guarantee that the same failure will not 
occur in the higher-order truncation. On the other hand, a term-by-term 
computation of the whole set of coefficients a n and b n is a quite nontrivial 
task due to a nonlinearity in the equations (28). 

Fortunately, another route is possible. Multiplying both the equations 
in (28) with (— k 2 ) n+1 , and performing a formal summation in n from zero 
to infinity, we arrive at the following expressions: 

{r oo n \ 

+ E E a n . m (-k 2 ) n - m a m (-k 2 ) m , (30) 
n=0m=0 J 

oo n 

B-b = A-a -k 2 Y / T, a n-m(-k 2 ) n - m b m (-k 2 r. 

n=0 m=0 

Now we notice that 

N n 

Jim Yl E a n „ m (-k 2 ) n - m a m (-k 2 r = A 2 , (31) 

N n 

lim E a n . m (-k 2 ) n - m b m (-k 2 r = AB. 

N ^°° n — n 
n=0 m=U 
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Accounting the initial condition (29), we come in (30) to a pair of coupled 
quadratic equations for the functions A and B: 

A = -\-k^-B + A*y (32) 

B = A{l-k 2 B). 

The result (32) concludes essentially the question of the computation of 
the functions A and B (22). Still, further simplifications are possible. In 
particular, it is convenient to reduce the consideration to a single function. 
Resolving the system (32) with respect to the function B, and introducing 
a new function, X(k 2 ) = k 2 B(k 2 ), we come to an equivalent cubic equation: 

_§<*_!).(*+«) = * (33) 

Since the functions A and B (22) are real-valued, we are interested only in 
the real- valued roots of eq. (33). 

An elementary analysis of this equation brings the following result: the 
real-valued root X(k 2 ) of the equation (33) is unique and negative for all 
finite values of the parameter k 2 . Moreover, the function X(k 2 ) is a mono- 
tonic function of k 2 . The limiting values are: 

lim X(k 2 ) = 0, lim X(k 2 ) = -0.8. (34) 

|fe|— >0 |fe|^oo 

The function X(k 2 ) is plotted in Fig. 1. 

Under the circumstances just mentioned, a function under the root in 
the Eq. (23) is negative for all values of the wave vector k, including the 
limits, and we come to the following dispersion law: 




, 5A2-16A + 20 
" ±= 2(T^A) ±Z tV 3 ' (35) 

where X = X(k 2 ) is the real- valued root of the equation (33). Since the func- 
tion X(k 2 ) is a negative function for all \k\ > 0, the attenuation rate, Reu±, 
is negative for all \k\ > 0, and the exact acoustic spectrum of the Chapman- 
Enskog procedure is stable for arbitrary wave lengths. In the short-wave 
limit, expression (35) reads: 

lim u± = -- ±i\k\V3. (36) 

fe|^oo 9 
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The characteristic equation of the original Grad equations (10) reads: 

3lo 3 + 3lo 2 + 9k 2 u + 5k 2 = 0. (37) 

The two complex-conjugated roots of this equation correspond to the hy- 
drodynamic modes, while the non-hydrodynamic real mode, u n h(k), and 
^nh(O) = —1, andw„/, — ► —0.5, as |fc| — > oo. Recall that the non-hydrodynamic 
modes of the Grad equations are characterized by the common property 
that for them u(0) / 0. These modes are irrelevant to the Chapman-Enskog 
method. As the final comment here, the expression (36) demonstrates a ten- 
dency of the exact attenuation rate, Re u>±, to a finite value, — | ~ —0.22, as 
\k\ — > oo. This asymptotics is in a complete agreement with the data for the 
hydrodynamic branch of the spectrum (37) of the original Grad equations 
(10). The attenuation rates (real parts of the dispersion relations u± for the 
Burnett (18), the super-Burnett (19), the exact Chapman-Enskog solution 
(35), are compared to each other in Fig. 2. In this figure, we also represent 
the attenuation rates of the hydrodynamic and of the non-hydrodynamic 
mode of the Grad equations (37). The results of this section lead to the 
following discussion: 

(i) The example considered above gives an opportunity to treat the prob- 
lem of a summation of the Chapman-Enskog expansion. The exact disper- 
sion relation (35) of the Chapman-Enskog procedure is demonstrated to be 
stable for all wave lengths, while the Bobylev instability is present on the 
level of the super-Burnett approximation. Moreover, it can be demonstrated 
that the function X (the real root of the equation (33)) is a real- valued ana- 
lytic function of the variable k. Thus, the treatment of the formal expansions 
performed above is justified. 

(ii) The exact result of the Chapman-Enskog procedure has a clear non- 
polynomial character. Indeed, this follows directly from (34): the function 
X(k 2 ) cannot be a polynomial because it maps the axis k into a segment 
[0, —0.8]. As a conjecture here, the resulting exact hydrodynamics is essen- 
tially nonlocal in space. For this reason, even if the hydrodynamic equa- 
tions of a certain level of the approximation is stable, it cannot reproduce 
the non-polynomial behavior for sufficiently short waves. 

(iii) The result of this section demonstrates that, at least in some cases, 
the sum of the Chapman-Enskog series amounts to a quite regular function, 
and the " smallness" of the Knudsen number e used to develop the Chapman- 
Enskog procedure (12) is no more necessary at the outcome. 
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2.2 The 3D10M Grad equations 

In this section we generalize our considerations of the Chapman-Enskog 
method the three-dimensional linearized 10-moment Grad equations [8]. The 
Chapman-Enskog series for the stress tensor, which is again due to a non- 
linear procedure, will be summed up in a closed form. The method we use 
follows essentially the one discussed above, though the computations are 
slightly more extensive. The reason to consider this example is that we 
would like to know what happens to the diffusive hydrodynamic mode in 
the short-wave domain. 

Throughout the section, we use the variables (1), and p and u are di- 
mensionless deviations of the pressure and of the mean flux from their equi- 
librium values, respectively. The point of departure is the set of the three- 
dimensional linearized Grad equations for the variables p, u, and cr, where 
cr is a dimensionless stress tensor: 

d t p = -|v-«, (38) 
dtu = —Vp — V ■ cr, 

dtcr = —Vu — -cr. 

e 

Equation (38) provides a simple model of a coupling of the hydrodynamic 
variables, u and to the non- hydrodynamic variable cr. These equations are 
suitable for an application of the Chapman-Enskog procedure. Therefore, 
our goal here is not to investigate the properties of eq. (38) as they are, but 
to reduce the description, and to get a closed set of equations with respect 
to the variables p and u alone. That is, we have to express cr in terms of 
spatial derivatives of p and of u. The Chapman-Enskog method, as applied 
to eq. (38) results in the following: 

oo 

cr = ]Te n+ V n ). (39) 

?1=0 

The coefficients cr*™) are due to the following recurrence procedure: 

n— 1 

ff W = - cfV™- 1 -™), (40) 

m=0 

where the Chapman-Enskog operators d[ m ^ act on the functions p and u, 
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and on their derivatives, as follows: 



«(™)n J ~ D ^P, m = 

a * DU = -DV- ff M), m >i ' W 



M = | -fDV-u, m = 
1 10, m > 1 

Here D is an arbitrary differential operator D = 11?= l $i > while ij is an 
arbitrary integer, and = 1. Finally, <x(°) = -Vu, which leads to the 
Navier-Stokes approximation. 

Our goal is to sum up the series (39) in a closed form. Firstly, we will 
make some preparations. 

As was demonstrated in [11], cr^ in the equations (39), (40), and (41), 
have the following explicit structure for arbitrary order n > (a generaliza- 
tion of the expressions (20) onto the three-dimensional case): 

cr( 2 ™) = o n A n V« + 6 n A n - 1 GV-«, (42) 
cr( 2 " +1 ) = c n A n G P , 

where A = V • V is the Laplace operator, and the operator G has the form: 

G = VV - -I A = -W. (43) 

The real- valued and yet unknown coefficients a n , b n , and c n in the equa- 
tion (42) are due to the recurrence procedure (40), and (41). Knowing the 
structure of the coefficients of the Chapman-Enskog series (42) , we can refor- 
mulate the Chapman-Enskog solution in terms of a self-consistent recurrence 
procedure for the coefficients a n , b n , and c n . Let us consider this derivation 
in more detail. 

The point of departure is the Fourier representation of the recurrence 
equations (40), (41), and (42). Writing 

u = U/, exp(ifc • x), 
p = p k exp(ik ■ x), 

and introducing the unity vector e k directed along k (k = ke k ), we come in 
the equations (40), (41), and (42) to the following: 

n-l 

(n) _ \ -> (n-l-m) 



- e ot^r m \ (44) 



m=0 
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flMn J -D k ikp k , m = 

9t DkUk = \ -D&.Jr*, ™>1 > (45) 
= ^-lD k ik.u k , m = 

where is an arbitrary tensor = Ils=i(*^s) is ! an d 

of n) = (-fe 2 ) n (a n ^ + 6 n i^(fc- M )), (46) 
= c n (—k 2 ) n+1 g k pk, 

where 

fffc = ( e fc e fc - = 2^^- ( 47 ) 

From the form of the Navier-Stokes approximation, a k °^ , it follows that 
ao = — 1 and 60 = 0, while a direct computation of the Burnett approxima- 
tion leads to: 

<r[ l) = \k 2 g kPk . (48) 

Thus, we have cq = —\, and with this the ansatz (42) is proven for n = 
in both the even and the odd orders. 

The further derivation relies upon induction. Let the structure (46) be 
proven up to the order n. The computation of the next, n+1 order coefficient 
<T fc 2 ^ n+1 ^ ' involves only the terms of the lower order. From the equation (44) 
we obtain: 

CT (2(n + l)) = ^(O^n+l) _ 2 ^ d (m )(T (2n + l-m)_ (4g) 

m=l 

The first term in the right hand side depends linearly on the coefficients c n : 

5 

= ^c n (-k 2 ) n+1 ig k k ■ u k . 

The remaining terms on the right hand side of the equation (49) contribute 
nonlinearly. Splitting the even and the odd orders of the Chapman-Enskog 
operators d^ m \ we rewrite the sum in the equation (49): 

- d { t m) cr k 2n+1 ' m) = - d[ 2l) a[ 2{n ~ l)+l) - d {2l+l) a k 2{n ' l)) (51) 

m=l 1=1 1=0 
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Due to (46) and (45), each term in the first sum is equal to zero, and we are 
left only with the second sum. We compute: 



^a+D^-O) = { - k y-\a n ^ikdf l+1) u k + b n ^ig k k ■ df + ^u k ), (52) 
while 

d\ 21 +l) u k = -(-k 2 ) l+1 ( ai u k + ±(a, + 2b l )e k (e k ■ u k )). (53) 
In the last expression, the use of the following identities was made: 

kku^ = k 2 (u k + ^e k (e k -u k )), (54) 
2 

kg k = -k. 

Substituting the expression (53) into the right hand side of the equation 
(52), and thereafter substituting the result into the right hand side of the 
expression (51), we come to the following in the right hand side of the 
equation (49): 



r (2(n + l)) = (.^l^^j^^.^l^ 

^2 \\{2a n -m + b n ^ m )(a m + 2b m ) + a n - m b m X ) ig k (k 



(55) 



• u k ). 



The functional structure of the right hand side of this expression is the same 
as that of the first of the equations in the set (46) , and thus we come to the 
first recurrence equation: 

a n+ iku k + b n+ ig k {k ■ u k ) = ^ 

ku k + 

\m=0 / 



c n + |^(2on-m + b n - m ){a m + 2b m ) + a n _ m 6 m |^j g k (k 



Considering in the same way the coefficient (j^ n+1 ^ +1 ^ we come to the 
second recurrence equation, 

2 n 

Cn+i = 2a n+ i + b n+1 + - Y ( 2a + bn—rn)c m . (57) 

O „ 

rn=0 
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Thus, the complete set of the recurrence equations is given by Eq. (56) and 
Eq. (57). Equation (56) is equivalent to a pair of scalar equations. Indeed, 
introducing new variables, 



gCn, (58) 



2 

Qn = -(2a n + b n ), 



and using the identity, 



ku k = (ku k - 2g k (k ■ u k )) + 2g k (k ■ u k ), 
and also noticing that 



g k : (ku k - 2g k {k ■ u k )) = 0, 

where : denotes the double contraction of tensors, we arrive in the equations 
(56) and (57) at the following three scalar recurrence relations in terms the 
coefficients r n , q n , and a n : 

n 

r n+1 = q n+ i + 

(59) 



m=0 



Qn+1 



5 

TT^n + ^ ' Qn—mQn 
6 m=0 



O-rt+1 — ^ ' O'n—m^m 

The initial condition for this system is provided by the explicit form of the 
Navier-Stokes and of the Burnett approximations, and it reads: 

T-o = -4/3, g = -4/3, ao = -1. (60) 

The recurrence relations (59) are completely equivalent to the original 
Chapman-Enskog procedure (40) and (41). In the one-dimensional case, 
the recurrence system (59) reduces to the first two equations for r n and 
q n . In this case, the system of recurrence equations is identical (up to the 
notations) to the recurrence system (28), considered in the preceding section. 
For what follows, it is important to notice that the recurrence equation for 
the coefficients a n is decoupled from the equations for the coefficients r n and 

Qn- 
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Now we will express the Chapman-Enskog series of the stress tensor (39) 
in terms of the coefficients r n , q n , and a n . Using again the Fourier transform, 
and substituting the expression (42) into the right hand side of the equation 
(39), we derive: 

cr k = A(k 2 )(klT k - 2g k {k ■ u fc )) + ^Q(k 2 )g k (k ■ u k ) - h 2 R(k 2 )g kPk , (61) 

From here on, we use a new spatial scale which amounts to k' = ek, and 
drop the prime. The functions A(k 2 ), Q(k 2 ), and R{k 2 ) in the equation (61) 
are defined by the power series with the coefficients due to (59): 

oo 

A(k 2 ) = Y,an(-k 2 T, (62) 

n=0 

oo 

Q(k 2 ) = ^Qni-^r, 

n=0 

oo 

R(k 2 ) = ]Tr n (-fc 2 r. 

n=0 

Thus, the question of summation of the Chapman-Enskog series (39) 
amounts to finding the three functions, A(k 2 ), Q{k 2 ), and R(k 2 ) (62) in the 
three- and two-dimensional cases, or to the two functions, Q(k 2 ), and R(k 2 ) 
in the one-dimensional case. 

Now we will focus on a problem of a computation of the functions (62) 
from the recurrence equations (59). At this point, it is worthwhile to notice 
again that a truncation at a certain n is not successful. Indeed, already in 
the one-dimensional case, retaining the coefficients qo, and ro, and q\ gives 
the super-Burnett approximation (16) which has the short-wave instability 
for k 2 > 3, as it was demonstrated in the preceding section, and there is no 
guarantee that the same will not occur in a higher-order truncation. 

Fortunately, the route of computations introduced in the preceding sec- 
tion works again. Multiplying each of the equations in (62) with (— k 2 ) n+1 , 
and performing a summation in n from zero to infinity, we derive: 

{r oo n \ 

3^ + E E Qn- m (-k 2 ) n - m q m (-k 2 ) m , (63) 
n=0 m=0 ) 

oo n 

R-r = Q-qo~k 2 J2 E qn- m (-k 2 ) n - m r m (-k 2 r, 

n=0 m=0 

oo n 

A-a = -k 2 J2J2 a n~m(-k 2 ) n - m a m (-k 2 r. 

n=0 m=0 
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Now we notice that 



N n 

lim E a n ^ m (-k 2 ) n - m a m (-k 2 r = A 2 , (64) 

N— >CO n — n 
n=U m=0 

AT n 

lim £ E Qn-m(—k 2 ) n ~ m r m (—k 2 ) m = QR, 

n=U m=0 

AT n 

lim V V q„{-k 2 ) n - m q m {-k 2 ) m = Q 2 . 

n=U m=0 

Taking into account the initial conditions (60) , and also using the expressions 
(64) , we come in the equation (63) to the following three quadratic equations 
for the functions A, R, and Q: 



Q = ~--k 2 ^-R + Q 2 j, (65) 

R = Q(l-k 2 R), 
A = -(l + k 2 A 2 ). 

The result (65) concludes essentially the question of computation of func- 
tions (62) in a closed form. Still, further simplifications are possible. In 
particular, it is convenient to reduce a consideration to a single function 
within the first two equations in the system (65). Introducing a new func- 
tion, X(k 2 ) = k 2 R(k 2 ), we come to an equivalent cubic equation: 

-\(X-l) 2 (X + ^) = ^. (66) 

This equation coincides with the equation (66) of the previous section. We 
will also rewrite the third equation in the system (65) using a function 
Y{k 2 ) = k 2 A{k 2 ): 

Y(l + Y) = -k 2 . (67) 

The functions of our interest (62) can be straightforwardly expressed in 
terms of the relevant solutions to the equations (66) and (67). Since all 
the functions in (62) are real-valued functions, we are interested only in the 
real- valued roots of the algebraic equations (66) and (67). 

The relevant analysis of the cubic equation (66) was already done above: 
the real- valued root X(k 2 ) is unique and negative for all finite values of 
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parameter k 2 . Limiting values of the function X{k 2 ) at k —> and at 
k — ► oo are given by the expression (34): 

limX(k 2 )=0, limX(k 2 )=--. 

k— >0 fc^oo 5 

The quadratic equation (67) has no real- valued solutions for k 2 > \, and 
it has two real- valued solution for each k 2 , where k 2 < \. We denote k c = \ 
as the corresponding critical value of wave vector. For k = 0, one of these 
roots is equal to zero, while the other is equal to one. The asymptotics 
Y — > 0, as k — > 0, answers the question which of these two roots of eq. (67) 
is relevant to the Chapman-Enskog solution, and we derive: 



\ (l - VI - 4/fc 2 ) k < k c 



y= 2^ v. ~, p«c (68) 
^ none k > k c 

The function Y (68) is negative for k < k c . 

From now on, X and Y will denote the relevant roots of the equations 
(66) and (67) just discussed. The Fourier image of the expression V • a 
follows from (61): 

ik- a k = Y((e k • u k )e k - u k ) - - - — (e k • u k )e k - iXkp k . (69) 

The latter expression contributes to the right-hand side of the second of 
equations in the Grad system (38) (more specifically, it contributes to the 
corresponding Fourier transform of this equation). Knowing (69), we can 
calculate the dispersion u(k) of the plane waves ~ expjwi + ik ■ x} which 
now follows from the exact solution of the Chapman-Enskog procedure. The 
calculation of the dispersion relation amounts to an evaluation of the deter- 
minant of a (d + 1) x (d + 1) matrix, and is quite standard (see, e.g. [21]). 
We therefore reproduce only the final result. The exact dispersion relation 
of the hydrodynamic modes reads: 

\d-i ( , ,2 X , , , 5 7,2/ 



( W - Y) a - 1 \u 2 - -^u; + -k\l -X))=0. (70) 

Here d is the spatial dimension. 

From the dispersion relation (70), we easily derive the following classifi- 
cation of the hydrodynamic modes: 

(i) For d = 1, the spectrum of the hydrodynamic modes is purely acoustic 
with the dispersion cj a which is given by the expression (35): 
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, , 5X2-16X + 20 
^ = —1 V^ ±l n\l § > ( 71 ) 

where X = X(/c 2 ) is the real- valued root of eq.(66). Since X is a negative 
function for all k > 0, the attenuation rate of the acoustic modes, Re u a , 
is negative for all k > 0, and the exact acoustic spectrum of the Chapman- 
Enskog procedure is free of the Bobylev instability for arbitrary wave lengths. 

(ii) For d > 1, the dispersion of the acoustic modes is given by the equa- 
tion (71). As follows from the Chapman-Enskog procedure, the diffusion- like 
(real- valued) mode has the dispersion u&: 



-\ (l- Vl-4fc 2 ) k<k c 
none k > k c 



u, d = n x vx ™) (72) 



The diffusion mode is (d — 1) times degenerated, the corresponding attenu- 
ation rate is negative for k < k c , and this mode cannot be extended beyond 
the critical value k c = \ within the Chapman-Enskog method. 

The reason why this rather remarkable peculiarity of the Chapman- 
Enskog procedure occurs can be found upon a closer investigation of the 
spectrum of the underlying Grad moment system (38). 

Indeed, in the original system (38), besides the hydrodynamic modes, 
there exist several non-hydrodynamic modes which are irrelevant to the 
Chapman-Enskog solution. All these non-hydrodynamic modes are charac- 
terized with a property that corresponding dispersion relations uj(k) do not 
go to zero, as k — > 0. In the point k c = ^, the diffusion branch (72) inter- 
sects with one of the non-hydrodynamic branches of the equation (38). For 
the larger values of the wave vector k, these two branches produce a pair of 
complex conjugated solutions with the real part equal to — \. Thus, though 
the spectrum of the original equations (38) continues indeed after k c , the 
Chapman-Enskog method does not recognize this extension as a part of the 
hydrodynamic branch. It is also interesting to notice that if we would accept 
all the roots of the equation (67), including the complex values for k > k c , 
and not only the real-valued root as was suggested by the asymptotics of 
the Chapman-Enskog solution (see the explanations preceding the equation 
(68)), then we would come in (70) to the structure of the dispersion relation 
just mentioned. 

The attenuation rates (the functions Re u a and Re ujd) are represented 
in Fig. 3, together with the relevant dependencies for the approximations of 
the Chapman-Enskog method. The non-hydrodynamic branch of eq. (38) 
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which causes the breakdown of the Chapman-Enskog solution is also repre- 
sented in Fig. 3. It is rather remarkable that while the exact hydrodynamic 
description becomes inapplicable for the diffusion branch at k > k c , the 
usual Navier-Stokes description is still providing a good approximation to 
the acoustic mode around this point. 

The analysis of this section leads to the following additional remarks to 
the conclusions made in the end of the section 2.1.: 

(i) The example considered above gives an opportunity to discuss the 
features of Chapman-Enskog solutions and the problem of an extension of 
the hydrodynamic modes into a highly non-equilibrium domain on the exact 
basis and in the full spatial dimension. The exact acoustic mode in the 
framework of the Chapman-Enskog procedure is demonstrated to be stable 
for all wave lengths, while the diffusion-like mode can be regarded for the 
hydrodynamic mode only in a bounded domain k < k c . It is remarkable that 
the result of the Chapman-Enskog procedure has a clear non-polynomial 
character. As a conjecture here, the resulting hydrodynamics is essentially 
nonlocal in space. It is also clear that any polynomial approximation to the 
Chapman-Enskog series will fail to reproduce the peculiarity of the diffusion 
mode demonstrated in frames of the exact solution. 

(ii) Concerning an extension of hydrodynamics into a highly non-equilibrium 
domain on the basis of the Boltzmann equations, this question remains open 

in a sense of an exact summation as above. In this respect, results for simpli- 
fied models can serve either for a test of approximate procedures or at least 
for a guide. In particular, the mechanism of the singularity of the diffusion- 
like mode through a coupling to the non-hydrodynamic mode might be a 
rather general mechanism of limitation of the hydrodynamic description, 
and not just a feature of the Grad systems. 

(iii) The result of this section demonstrates that the sum of the Chapman- 
Enskog series amounts to either a quite regular function (as is the function 
X), or to a function with a singularity at finite k c . In both cases, however, 
the "smallness" of the Knudsen number e used to develop the Chapman- 
Enskog procedure plays no role in the outcome of the Chapman-Enskog 
procedure. 
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3 The dynamic invariance principle 



3.1 Partial summation of the Chapman-Enskog expansion 

The examples considered above demonstrate that it makes sense to speak 
about the sum of the Chapman-Enskog expansion, at least when the Chapman- 
Enskog method is applied to the (linearized) Grad equations. However, even 
in this case, the possibility to perform the summation exactly seems to be 
a lucky exceptions rather than a rule. Indeed, computations become more 
bulky with the increase of the number of the moments included in the Grad 
equations. Therefore, we arrive at a question: how can we approximate the 
recurrence equations of the Chapman-Enskog method to account all the or- 
ders in the Knudsen number? Any such method amounts to some "partial" 
summation of the Chapman-Enskog expansion, and this type of working 
with formal series is widely spread in various fields of physics. 

In this section we will discuss a method of approximating the Chapman- 
Enskog expansion in a whole. As we now have the exact expressions for the 
Chapman-Enskog solution for the linearized 10 moment Grad equations, it 
is natural to start with this example for the reason of comparison. 

Let us come back to the originating one-dimensional Grad equations 
(10), and to the corresponding formulas of the Chapman-Enskog method 
(12) and (13). Instead of using the exact equations (12) in each order n, we 
introduce the following approximate equations: 

Let N > 1 is some fixed integer. Then, instead of equations (12), we 
write: 

n-l 

a (n) = ^M^n-i-^ n<N, (73) 

m=0 
N-l 

a (n) = _Y^d\ m) a { - n - 1 - m \ n>N. (74) 

m=0 

This approximation amounts to the following: up to the order N, the 
Chapman-Enskog procedure (12) is taken exactly (equation (73)), while in 
the computations of the higher orders (equation (74)) we restrict the set 
of the Chapman-Enskog operators (13) only to the order N. Thus, the 
Chapman-Enskog coefficients of the order higher than N are taken into 
account only "partially". As N tends to infinity, the recurrence procedure 
(73) and (74) tends formally to the exact Chapman-Enskog procedure (12). 
We will further refer to the equations (73) and (74) as to the regularization of 
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the N — th order. In particular, taking N = 1, we come to the regularization 
of the Burnett approximation, taking N = 2 we come to the regularization 
of the super-Burnett approximation, etc. 

It can be demonstrated that the approximate procedure just described 
does not alter the structure of the functions <j( 2ra ) and a^ 2n+1 ^ (20), while the 
recurrence equations for the coefficients a n and b n (20) will differ from the 
exact result of the full Chapman-Enskog procedure (28). The advantage of 
the regularization procedure (73) and (74) above the exact Chapman-Enskog 
recurrence procedure (12) is that the resulting equations for the coefficients 
a n and b n are always linear, as they appear from the equations (73) and 
(74). This feature enables one to sum up the corresponding series exactly, 
even if the originating nonlinear procedure leads to a too difficult analysis. 
The number TV can be called the "depth" of the approximation: the greater 
is N the more low-order terms of the Chapman-Enskog expansion are taken 
into account exactly due to the equations (73). 

For the first example, let us take N = 1 in the equations (73) and 
(74). The regularization of the Burnett approximation in accord with these 
equations reads: 

a (n) = _ d (0) a (n-l)^ (75) 

where n > 1, and er(°) = — (4/3)d x ii. Turning to the Fourier variables, we 
derive: 

af n) = a n (-k 2 ) n iku k , (76) 

4 2n+1) = K(-k 2 r +l Pk , 

while the coefficients a n and b n are due to the following recurrence procedure: 
a-n+i = -^K, b n = a n , ao = —-, (77) 

and whereupon 

a n = K = ( - J a . (78) 

Thus, denoting as a^ k the Fourier transform of the regularized Burnett ap- 
proximation, we obtain: 

a ^ = - 3+ 4 5fc2 { iku k ~ k2 Pk) ■ (79) 

It should be noticed that the recurrence equations (77) can also be ob- 
tained from the exact recurrence equations (28) by canceling the nonlinear 
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terms. Thus, the approximation adopted within the regularization proce- 
dure (75) amounts to the following rational approximation of the functions 
A and B (22): 

i4f - Bf — 5W <80) 

Substituting the latter expressions instead of the functions A and B in the 
formula for the dispersion (23), we come to the dispersion relation of the 
hydrodynamic modes within the regularized Burnett approximation: 



2k 2 ,,, 75k 2 k 2 + 66k 2 + 15 . , 

W± = -3T^ ±Akl i 25fe2fc2 + 3 0fc 2 + 9 - ( 81 ) 

The dispersion relation (81) is stable for all wave vectors, and in the short- 
wave limit we have: 

lim uj ± = -0.4±i|fc|V3. (82) 

Thus, the regularized Burnett approximation leads to a qualitatively the 
same behavior of the dispersion relation, as it is for the exact result (36), 
giving the limiting value of the attenuation rate equal to —0.4 instead of the 
exact value —2/9. 

Consider now the regularization of the super-Burnett approximation. 
This amounts to setting N = 2 in the recurrence equations (73) and (74). 
Then, instead of the equations (75), we have: 

a (D = _ flt (o) ff (o )> (g3) 

a (2+n) = _ a (0) <T (n+l)_ a (l) <T (n) > 

where n > 0. The corresponding recurrence equations for the coefficients a n 
and b n now are as follows: 

1 4 

a n +i = g&n, a n = b n , a = --. (84) 

Thus, instead of the expressions (80), we come to the following: 

A * ~ B * ~ -5TF- <85) 

The corresponding dispersion relation of the regularized super-Burnett ap- 
proximation reads: 



2k 2 , ,,. 25k 2 k 2 + 78k 2 + 45 

w± = ~sT¥ ± ,|fe| VWTl8M' (86) 
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while in the short-wave limit the following asymptotics takes place: 



The Bobylev instability is removed again within the regularization of 
the super-Burnett approximation, and also the lower-order terms of the 
Chapman-Enskog expansion are taken into account more precisely in com- 
parison to the regularized Burnett approximation. However, the approxima- 
tion in a whole has not improved (see Fig. 4). Thus, we can conclude that 
though the partial summation method (73) and (74) is in a capacity to remove 
the Bobylev instability, and to reproduce qualitatively the exact Chapman- 
Enskog solution in the short-wave domain, the exactness, however, does not 
increase monotonically with the depth of the approximation N. This draw- 
back of the regularization procedure indicates once again that an attempt to 
capture the lower-order terms of the Chapman-Enskog procedure does not 
succeed in a better approximation in a whole. 

3.2 The dynamic invar iance 

All the procedures considered so far (exact or approximate) were taking the 
Chapman-Enskog expansion as the starting point. However, the result of 
the summation in these procedures does not involve the Knudsen number 
e explicitly, neither the sum does apply for a "smallness" of this parame- 
ter. Therefore, it makes sense to reformulate the problem of the reduced 
description (for the Grad equations (10) this amounts to the problem of 
constructing a function ak(uk,Pk, k)) in a way where the parameter e does 
not come into play at all. Further, in a framework of such an approach, we 
can seek for a method of explicit constructing the function o- k (u k ,Pk, k), and 
which does not rely upon the Taylor-like expansions as above. 

In this section we introduce the approach just mentioned, considering 
again the illustrative example (10). These ideas will be extensively used 
in the sequel, and they also constitute the basis of the so-called method of 
invariant manifold for dissipative systems [23]. 

Let us rewrite here the equations (10), in the Fourier variables, and 
canceling the parameter e: 




(87) 



dtPk = ~^iku k , 
d t u k = -ikp k - ika k , 



(88) 
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d t a k = ~^iku k - o-fe. 

The result of the reduction in the system (88) amounts to a function 
a k (u k ,p k , k), which depend parametrically on the hydrodynamic variables 
u k and pk, and also on the wave vector k. Due to the linearity of the problem 
under the consideration, this function depends linearly on u k and p k , and 
we can start with the form given by the equation (21): 

Vk(uk,Pk, k) = ikAu k - k 2 Bp k , (89) 

where A and B are undetermined functions of k. Now, however, we do not 
refer to a power series representation of these functions as in the equations 
(22). 

Given the form of the function cr k (u k ,p k , k) (89), we can compute its 
time derivative in two different ways. On the one hand, substituting the 
function (89) into the right hand side of the third of the equations in the set 
(88), we derive: 

df CTO a k = -ik (| + A) u k + k 2 Bp k . (90) 

On the other hand, computing the time derivative due to the first two equa- 
tions (88), we come to the following: 

dr cro ak = plQ tUk + pl dtpk (91) 

du k dp k 



= ikA (—ikp k — iko~ k ) — k B ^— -iku k 

= ik {^k 2 B + k 2 A^j u k + k 2 ( K A- k 2 B) p k . 

Equating the expressions in the right hand sides of the equations (90) and 
(91), and requiring that the resulting equality holds for any values of the 
variables u k and p k , we come to the following two algebraic equations: 

F(A, B, k) = -A-\-k 2 ( b -B + A 2 ]=^ (92) 



3 V3 

G(A, B, k) = -B + A(l-k 2 B) = 0. 

These equations are nothing else but the equations (32). Recall that equa- 
tions (32) have been obtained upon the summation of the Chapman-Enskog 
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expansion. Now, however, we have come to the same result without using 
the expansion. Thus, the equations (92) (or, equivalently, (32)) can be used 
as a starting point for the constructing the function (89). 

Now it is important to comment on the somewhat formal manipulations 
which have led to the equations (92). First of all, by the very sense of the 
reduced description problem, we are looking for a set of functions which 
depend on the time only through the time dependence of the hydrodynamic 
variables and p^- That is, we are looking for a set (89), which is pa- 
rameterized with the values of the hydrodynamic variables. Further, the 
two time derivatives, (90) and (91), are relevant to the "microscopic" and 
the "macroscopic" evolution within the set (89), respectively. Indeed, the 
expression in the right hand side of the equation (90) is just the value of the 
vector field of the original Grad equation in the points of the set (89). On 
the other hand, the expression (91) reflects the time derivative due to the 
reduced (macroscopic) dynamics, which, in turn, is self-consistently defined 
by the form (89). The equations (92) provide, therefore, the dynamic in- 
variance condition of the reduced description for the set (89): the function 
o~k{uk{t),Pk{t), k) is a solution to both the full Grad system (88) and to the 
reduced system which consists of the two first (hydrodynamic) equations. 
For this reason, the equations (92) and their analogs which will be obtained 
on the similar reasoning, will be called the invariance equations. 

3.3 The Newton method 

Let us concentrate on the problem of solving the invariance equations (92). 
Clearly, if we are going to develop the functions A and B into the power 
series (22), we will come back to the Chapman-Enskog procedure. Now, 
however, we see that the Chapman-Enskog expansion is just a method to 
solve the invariance equations (92), and maybe not the optimal one. 

Another possibility is to use iterative methods. Indeed, let us apply 
the Newton method. The algorithm of is as follows: Let A$ and Bq are 
some initial approximations chosen for the procedure. The correction, A\ = 
Aq + SAi and B\ = Bq + 5Bi, due to the Newton iteration is obtained upon 
a linearization of the equations (92) about the approximation Aq and Bo. 
Computing the derivatives, we can represent the equation of the Newton 
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iteration in the matrix form: 

/ dF(A,B,k) | dF(A,B,k) 



dkj \a=Ao,b=b mz~ \a=a ,b=b I f F ( A o,B ,k) 

(93) 



9 -^A m \A^B ^^L=. ,S=S M Tl G(4,,*>,*) 




where 



a=a ,b=b = -(l + 2k 2 A ), (94) 



dF(A,B,k) 
OA 

dF(A, B, k) 5 2 



^ U=Ao,B=_B 

dG(A,B,k) 



dB 



\A=A ,B= 



1 - A; 2 £ , 

- (l + A; 2 A)) • 



Solving the system of linear algebraic equations, we come to the first cor- 
rection 8A\ and 5B\. Further corrections are found iteratively: 

A n+1 = A n + 5A n+1 , (95) 
Bn+i = B n + SB n+ i, 

where n > 0, and 



-{l + 2k 2 A n ) -Ik 2 \ ( 5A n+1 \ F(A n ,B n ,k) 
1 - k 2 B n -(1 + k 2 A n ) ) \ 8B n+1 ) + I G(A n , B n , k) 



(96) 

Within the algorithm just presented, we come to the problem how to 
choose the initial approximation Aq and Bq. Indeed, the method (95) and 
(96) is applicable formally to any initial approximation, however, the con- 
vergence (if at all) might be sensitive to the choice. 

For the first experiment let us take the Navier-Stokes approximation of 
the functions A and B: 

A = B = ~ 

The outcome of the first two Newton iterations (the attenuation rates as 
they follow from the first and the second Newton iteration) are presented 
in Fig. 5. It is clearly seen that the Newton iterations converge rapidly to 
the exact solution for moderate k, however, the asymptotic behavior in the 
short-wave domain is not improved. 



= 0. 
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Another possibility is to take the result of the regularization procedure 
as presented above. Let the regularized Burnett approximation (80) is taken 
for the initial approximation, that is: 

A ° = A ? = -JT^- B » = B * = (97 » 

Substituting the expression (97) into the equation (95) and (96) for n = 
0, and after some algebra, we come to the following first correction: 

4(27 + 63A; 2 + 153fc 2 fc 2 + 125fc 2 fc 2 fc 2 ) 
1 ~ 3(3 + 5£: 2 )(9 + 9£: 2 + 67£; 2 A; 2 + 75£: 2 fc 2 fc 2 )' ( ' 

4(9 + 33A: 2 + 115fc 2 fc 2 + 75k 2 k 2 k 2 ) 
Bl ~ ~ (3 + 5£: 2 )(9 + 9fc 2 + 67£: 2 £: 2 + 75k 2 k 2 k 2 ) 

The functions (98) are not yet the exact solution to the equations (92) (that 
is, the functions F(A\, B±, k) and G(Ai , B\ , k) are not equal to zero for all 
k). However, substituting the functions A\ and B\ instead of A and B into 
the dispersion relation (23), we derive in the short-wave limit: 

lim w± = -- ±i\k\y/S. (99) 

fe|^oo 9 

That is, already the first Newton iteration, as applied to the regularized 
Burnett approximation, leads to the exact expression in the short-wave do- 
main. Since the first Newton iteration appears to be asymptotically exact, 
the next iterations improve the solution only for the intermediate values of 
k, whereas the asymptotic behaviour remains exact in all the iterations. The 
attenuation rates for the first and for the second Newton iterations with the 
initial approximation (97) are represented in Fig. 6. The agreement with 
the exact solution is excellent. 

A one more test is to take the result of the super-Burnett approximation 
(85) for the initial condition in the Newton procedure (96). As we know, 
the regularization of the super-Burnett approximation is provides a poorer 
approximation in comparison to the approximation (97), in particular, in 
the short-wave domain. Nevertheless, the Newton iterations do converge 
though less rapidly (see Fig. 7). 

The examples considered so far demonstrate that the Newton method, 
as applied to the invariance equations (92) is a more powerful tool in com- 
parison to the Chapman-Enskog procedure. It is also important that the 
initial approximation should be "properly chosen", and should reproduce 
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the features of the solution in a whole (not only in the long- wave limit), at 
least qualitatively. 

The best among the initial approximations considered so far is the reg- 
ularized Burnett approximation (97). We have already commented on the 
relation of this approximation to the invariance equations, as well as on its 
relation to the Chapman-Enskog procedure. The further important obser- 
vation is as follows: 

Let us choose the Euler approximation for the functions A and B, that 

is: 

A) = B = (100) 



The equation of the first Newton iteration (96) is very simple: 




(101) 



and 

• 4i = Bi = -3tW^ (102) 

Thus, the regularized Burnett approximation is at the same time the first 
Newton correction as applied to the Euler initial approximation. This prop- 
ertie distinguishes the regularization of the Burnett approximation among 
other regularizations. Now the functions (98) can be regarded as the second 
Newton correction as applied to the Euler initial approximation (100). 

Finally, let us examine the question what does the Newton method do 
in a case of singularities. As we have demonstrated in the previous section, 
the singularity of the diffusion-like mode occurs when this mode couples to 
a non-hydrodynamic mode of the 10 moment Grad system if the spatial 
dimension is greater that one. 

We leave it here without a proof that the invariance equation method as 
applied to the 10 moment Grad system (38) leads to the system of equations 
(65). We have already demonstrated what is the outcome of the Newton 
method as applied to the first two equations of this system (responsible for 
the acoustic mode and containing no singularities). The Newton method, 
as applied to the equation (67), reads: 

Y n+1 = Y n + SY n+1 , (103) 
{l + 2Y n )5Y n+1 + {Y n {l+Y n ) + k 2 } = 0, 
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where n > 0, and Yq is a chosen initial approximation. Taking the Euler 
approximation (Yq = 0), we derive: 

Fi = -k 2 , (104) 
fc 2 (l + fc 2 ) 
2 l-2k 2 ' 

The second approximation, Y^, is singular at the point &2 = \/l/2, and it 
can be demonstrated that all the further corrections do also have the first 
singularity at points k n , while the sequence k2, ■ ■ ■ ,k n tends to the actual 
branching point of the invariance equation (67) k c = 1/2. The analysis of 
further corrections demonstrates that the convergence is very rapid (see Fig. 
8). 

The approximations (104) demonstrate that unlike the polynomial ap- 
proximations, the Newton method is sensitive to detect the actual singu- 
larities of the hydrodynamic spectrum. Formally, the function Y2 becomes 
positive as k becomes larger than &2, and thus attenuation rate, uj^ = Y2 
becomes positive after this point. However, unlike the super-Burnett ap- 
proximation for the acoustic mode, this transition occurs now in a singular 
point. Indeed, the attenuation rate Y2 tends to "minus infinity", as k tends 
to ki from the left. Thus, as described with the Newton procedure, the 
non-physical domain is separated from the physical one with an "infinitely 
viscid" threshold. The occurrence of the poles in the Newton iterations is, 
of course, quite clear. Indeed, the Newton method involves the derivative 
of the function R(Y) = Y(Y + 1) which appears on the left hand side of 
the equation (67). The derivative dR(Y)/dY becomes zero in the singular- 
ity point Y c = —1/2. The results of this section bring us to the following 
discussion: 

(i) The result of the exact summation of the Chapman-Enskog procedure 
brings us to the same system of equations as the principle of the dynamic 
invariance. This is demonstrated above for a specific situation but it holds 
as well for any (linearized) Grad system. The resulting equations are always 
nonlinear (even for the simplest linearized kinetic systems, such as Grad 
equations). 

(ii) Now we are able to alter the viewpoint: the set of the invariance 
equations can be considered as the basic in the theory, while the Chapman- 
Enskog method is a way to solve it via an expansion in powers of k. The 
method of power series expansion is neither the only method to solve equa- 
tions, nor the optimal. Alternative iteration methods might be better suited 
to the problem of constructing the reduced description. 
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(iii) An opportunity to derive the invariance equation in a closed form, 
and next to solve it this or that way, is, of course, rather exotic. The situa- 
tion becomes complicated already for the nonlinear Grad equations, and we 
should not expect anything simple in the case of the Boltzmann equation. 
Therefore, if we are willing to proceed along these lines in other problems, 
the attention draws towards the approximate procedures. With this, the 
question appears: what amount of information is required to execute the 
procedures? Indeed, the Navier-Stokes approximation can be obtained with- 
out any knowledge of the whole nonlinear system of invariance equations. It 
is important that the Newton method, as applied to our problem, does not 
require any global information as well. This was demonstrated above by a 
relation between the first iteration as applied to the Euler approximation 
and the regularization of the Burnett approximation. 

3.4 Invariance equation for the 1D13M Grad system 

Let us consider as the next example the problem of the reduced description 
for the one-dimensional thirteen moment Grad system. Using the dimen- 
sionless variables as above, we write the one-dimensional version of the Grad 
equations (2) and (3) in the /c-representation: 



dtPk = 


-iku k , 




t u k = 


-ikpk - 


ikT k — ikok, 


d t T k = 


--iku k 


- ^ikqk, 


d t <Jk = 


-^iku k 


- T=ikqk ~ (Tk 
15 


d t qk = 


-\ikT k 


- ikuk - -qk- 



The Grad system (105) provides the simplest coupling of the hydrodynamic 
variables p k , Uk, and T k to the non- hydrodynamic variables, a k and qk, the 
latter corresponds to the heat flux. As above, our goal is to reduce the 
description for the Grad system (105) to the three hydrodynamic equations 
with respect to the variables p k , Uk, and T k . That is, we have to express the 
functions Ok and qk in terms of pk, Uk, and T^: 

Ofc = a k (pk,u k ,T k ,k), 
qk = qk(pk,Uk,T k ,k). 
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The Chapman-Enskog method, as applied for this purpose, results in the 
following algebraic scheme (we omit the Knudsen number e): 



' n-l 



(n) I \— * o( m ) (n—l—m) , « ., (n— 1) I /-, nn \ 

- \ 2^ d t a k + T7M k \ (106) 



(n) 



,m=0 

■ n-l 



15 



- -\Y. d< t m) it 1 - m) + ^t 1) 



km=0 J 

where the Chapman-Enskog operators act as follows: 

^ Pfc " \ 0, m > 1 ' (107) 

Q (m) f -*fc(Pk + T fc) "1 = 

* Ufc = \ m>l ' 



m = 

-^ikq^™ 1 ^, m > 1 

The initial condition for the recurrence procedure (106) reads: erj^ = —7-ikuk, 

and ^°' ) = — ^ifcTfe, which leads to the Navier-Stokes- Fourier hydrodynamic 
equations. 

Computing the coefficients and gj^, we come to the Burnett ap- 
proximation: 

(Tife = -^i/cu fc + ^k 2 p k - ^k 2 T k , (108) 

15 lm 7 l2 
9iJfc = -^-ikT k + -k 

The Burnett approximation (108) coincides with that obtained from the 
Boltzmann equation, and it is precisely the case where the instability was 
first demonstrated in the paper [2] . 

The structure of the terms and qj^ (an analog of the equations (20) 
and (42)) is as follows: 

4 2n) = a n (-k 2 ) n iku k , (109) 

4 2n+1) = b n (-k 2 r^ Pn+ c n (-k 2 r^T k , 

Qk n) = Pn(-k 2 ) n ikp k + 7n(-k 2 ) n ikiTk, 

qt +1) = M-k 2 r +1 uk. 
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A derivation of the invariance equation for the system (105) goes along 
the same lines as in the previous section. We seek the functions of the 
reduced description in the form: 

a k = ikAu k - k 2 B Pk - k 2 CT k , (110) 
q k = ikXp k + ikYT k - k 2 Zu k , 

where the functions A, . . . Z are a subject of a further analysis. 

The invariance condition results in a closed system of equations for the 
functions A, B, C, X, Y, and Z. As above, computing the microscopic time 
derivative of the functions (110), due to the two last equations of the Grad 
system (105) we derive: 

df CIO <r k = -ik[^-^k 2 Z + A^u k (111) 

+k 2 (J^X + B^jp k + k 2 (J^Y + Cj T k , 
8f cro q k = k 2 (A + -z)u k +ik(k 2 B --x) Pk -ik(- -k 2 C --yW. 



On the other hand, computing the macroscopic time derivative due to the 
first three equations of the system (105), we obtain: 

amacro ® a k ^ a kr, . d(J k 

d t <J k = -7T- d t u k + — d t p + — — d t T k (112 
ou k dp k dT k 

= ik(k 2 A 2 + k 2 B + hi 2 C -hi 2 k 2 CZ^u k 
+ (k 2 A - k 2 k 2 AB - \ 2 k 2 CX^ Pk 

+ (k 2 A-k 2 k 2 AC-h 2 k 2 CY^)T k ; 

omacro 9q k dq k dq k 

du k dp k dT k 

-k 2 k 2 ZA + k 2 X + \k 2 Y - \k 2 k 2 Yz\ u k 

o o J 



+ik [k 2 Z - k 2 k 2 ZB + ^k 2 YX 



Pk 



+ik (k 2 Z - k 2 k 2 ZC + \k 2 Y 2 ) 7 / 



39 



Equating the corresponding expressions in the formulas (111) and (112), we 
come to the following system of coupled equations: 



Fi = 


-t + ±k 2 Z -A- k 2 A 2 - k 2 B - -k 2 C + -k 2 k 2 CZ 
3 15 3 3 


F 2 = 


—X + B-A + k 2 AB + -k 2 CX = 0, 
15 3 


F 3 = 


— Y + C-A + k 2 AC + -k 2 CY = 0, 
15 3 


F 4 = 


A + -Z + k 2 ZA -X--Y + -k 2 YZ = 0, 
3 3 3 


Fs = 


k 2 B --X- k 2 Z + k 2 k 2 ZB - \k 2 YX = 0, 


F 6 = 


-| + k2 ° ~\ Y - k2z + k 2 k 2 ZC - \k 2 Y 2 = 0. 



As above, the invariance equations (113) can be also obtained upon the 
summation of the Chapman-Enskog expansion, after the Chapman-Enskog 
procedure is casted into a recurrence relations for the coefficients a n , . . . , a n 
(109). This route is less straightforward than the one just presented, and 
we omit the proof. 

The Newton method, as applied to the system (113), results in the fol- 
lowing algorithm: 

Denote as A the six-component vector function A = (A,B,C,X,Y,Z). 
Let Aq is the initial approximation, then: 

A n+1 = A n + SA n+1 , (114) 

where n > 0, and the vector function 5A n+ i is a solution to the linear system 
of equations: 

N n 6A n+1 + F n = 0. (115) 

Here F n is the vector function with the components Fi(A n ), and N n is a 
6x6 matrix: 



-{l + 2k 2 A n ) 
k 2 B n - 1 
k C n — 1 
1 + k 2 Z n 







-k 2 
1 + k 2 



k 2 {l + k 2 Z n ) 




-2/3fc 2 (l- k 2 Z n ) 

2/3k 2 X n 
l + 2/3k 2 Y n + k 2 A r , 



k 2 (l + k 2 Z n ) 



(116) 
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2/3k 2 (4/5 + k 2 C n ) 

2/3(4/5 + k 2 C n ) 

2/3(4/5 + k 2 C n ) 

-1 -2/3(1 - k 2 Z n ) 2/3 + k 2 A n + 2/3k 2 Y r 

-2/3(1 + k 2 Y n ) -2/3k 2 X n -k 2 {l - k 2 B n ) 

-2/3(1 + 2k 2 Y n ) -k 2 {l-k 2 C n ) ) 

The Euler approximation gives: Aq = . . . = Zq = 0, while Fi = —4/3, 
F 6 = -5/2, and F 2 = . . . = F 5 = 0. The first Newton iteration (115) 
as applied to this initial approximation, leads again to a simple algebraic 
problem, and we have finally obtained: 

141fc 2 + 20 

1 " ° 867F+2105F+300' ( ' 

459k 2 k 2 + 810A: 2 + 100 

B\ = —20 



3468k 2 k 2 k 2 + 12755/e 2 /c 2 + 11725/c 2 + 1500' 
51k 2 k 2 - 485A; 2 - 100 

Gi = — 10-; 



Zi = -15 



3468fc 2 fc 2 A; 2 + 12755fc 2 fe 2 + 11725& 2 + 1500' 

375fc 2 (21fc 2 - 5) 

2(3468/fc 2 A; 2 A; 2 + 12755/c 2 /c 2 + 11725A; 2 + 1500) ' 

225(394fc 2 /c 2 + 685A; 2 + 100) 
"4(3468A; 2 A; 2 fc 2 + 12755k 2 k 2 + 11725A: 2 + 1500) ' 
153fc 2 + 35 



867A; 4 + 2105A; 2 + 300 ' 



Substituting the expression (109) into the first three equations of the 
Grad system (105), and proceeding to the dispersion relation as above, we 
derive the latter in terms of the functions A,...,Z: 

^ - k 2 {^Y + A^u 2 (118) 

+ k 2 (- - -k 2 Z - -k 2 C - k 2 B + -k 2 AY + -k 2 k 2 Cz) to 
V3 3 3 3 3 J 

+ ^k 2 (k 2 X-k 2 Y + k 2 k 2 BY-k 2 k 2 XC) = 0. 

When the functions A±, . . . , Z\ (117) are substituted instead of A, . . . , Z 
into the equation (118), the dispersion relation of the first Newton iteration, 
as applied to the invariance equations (113) with the Euler initial approx- 
imation, is obtained. This result coincides with the regularization of the 
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Burnett approximation, which was considered in [11]. There it was demon- 
strate that the equilibrium is stable within this approximation for arbitrary 
wave lengths. The dispersion relation for the Burnett approximation, in 
turn, is due to the approximation 

A = -4/3, B = -4/3, C = 2/3, X = 0, Y = -15/4, Z = -7/4, 

as it follows from a comparison of the expressions (108) and (110). The 
dispersion relation for the Burnett approximation coincides with the one 
obtained in [2] from the Boltzmann equation. 



3.5 Invariance equation for the 3D13M Grad system 

The final example to be considered is the 13 moment Grad system in the 
full spatial dimension, (2) and (3). Let us rewrite here the original system 
in terms of Fourier variables: 

dtPk = -ike k ■ u k , (119) 

dtu k = -ike k p k - ike k T k - ike k ■ a k , 
2 

dtT k = --ik(e k -u k + e k -q k ), 

dt&k = -ike^Uk -^ike^ - a k , 

5 2 
dtQk = ~-ike k T k -ike k - a k - -q k . 

Here we have represented the wave vector k as k = ke k , and e k is the unit 
vector. 

The structure of the even and odd Chapman-Enskog coefficients, cr k n ^ 
and q k n \ turns out to be as follows: 

a k 2n) = (-k 2 ) n ik {a n (epi£ - 2g k (e k • u k )) + b n g k (e k • u k )} ,(120) 

Q { k n) = (-k 2 ) n ike k { ln T k + 5 nPk }, 
g (2n+i) = (-^^{a^ej^.^ . Ufe ) +^ n ( Ufc _ efc (e fc . Wfc )} ) 

where g k = l/2e k e k , and the real- valued coefficients a n ,...,/3 n are due to 
the Chapman-Enskog procedure (7) and (8). 
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The expressions just presented prompt that the dynamic invariant form 
of the stress tensor and of the heat flux reads: 



a k = ikA(e k u k -2g k (e k -u k )) + 2ikBg k (e k -u k ) (121) 

-2k 2 Cg k T k - 2k 2 Dg kPk , 
q k = ikZe k T k + ikUe k p k 

-k 2 X(u k - e k (e k ■ u k )) - k 2 Ye k (e k ■ u k ), 

where the functions A, . . . , Y depend on k. The dynamic invariance condi- 
tion results in the following two closed systems for these functions: 

-U + D- B+-k 2 CU +-k 2 BD = 0, (122) 
5 3 3 

-Z + C- B + -k 2 CZ + -k 2 BC = 0, 
5 3 3 

-1 + -k 2 Y - B - -k 2 C - k 2 D - -k 2 B 2 + -k 2 k 2 CY = 0, 
5 3 3 3 

i k i D _ lu _ k 2 Y - -k 2 ZU + -k 2 k 2 YD = 0, 
3 3 3 3 

J_ + i k z C --Z-k 2 Y--k 2 Z 2 + -k 2 k 2 YC = 0, 
2 3 3 3 3 

4 2 2 2 4 

-B + -Y - U - -Z + -k 2 ZY + -k 2 YB = 0, 

3 3 3 3 3 

and 

-l-A + -k 2 X-k 2 A 2 = 0, (123) 
5 

A + -X + k 2 AX = 
3 

The method of summation of the Chapman-Enskog expansion can be also 
developed, starting with the structure of the Chapman-Enskog coefficients 
(120), and in the same manner as in section 2. Simple but rather extensive 
computations in this case lead, of course, to the invariance equations (122) 
and (123). 

The Newton method, as applied to the systems (122) and (123) with the 
initial Euler approximation, leads in the first iteration to the regularization 
of the Burnett approximation reported earlier in [11]. 

Introducing the functions A = k 2 A and X = k 2 X, we come in the 
equation (123) to the following: 

V ; 4(6A + 5) ' V ' 
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while 



X = - 



3A 



2 + 3A' 



The derivative, dR(A)/dA, becomes equal to zero in the point A c rj —0.364, 



which gives the critical wave vector k c = y —R(A C ) ~ 0.305. The Newton 
method, as applied to the equation (124) with the initial Euler condition 
A = 0, gives the following: the outcomes of the first and of the second 
iterations are regular functions, while the third and the further iterations 
bring a singularity which converges to the point k c (see Fig. 9). These 
singularities (the real poles) of the Newton corrections are of the same nature 
as discussed above. 

3.6 Gradient expansions in kinetic theory of phonons 

3.6.1 Exact Chapman-Enskog solution and onset of second sound 

In this section, we close our discussion of linearized Grad systems with an 
application to simple models of the phonon transport in rigid insulators. 
It is demonstrated that the extended diffusion mode transforms into second 
sound mode due to its coupling to a non-hydrodynamic mode at some critical 
value of the wave vector. This criticality shows up as a branching point of 
the extension of the diffusion mode within the Chapman-Enskog method. 
Though the analysis is essentially similar to the examples considered above, 
it is presented in some details for the sake of completeness. 

Experiments on heat pulses propagation through crystalline media [24] 
confirmed existence of a temperature window (the Guyer-Krumhansl win- 
dow [25]) with respect to which the features of heat propagation are quali- 
tatively different: At temperatures exceeding the high-temperature edge of 
the window, the heat propagates in a diffusion-like way. Below the low- 
temperature edge of the window, the propagation goes in a ballistic way, 
with a constant speed of sound. Within the window, the propagation be- 
comes wave-like. This latter regime is called second sound (see [26] for a 
review) . 

This problem has drawn some renewed attention in the last years. Mod- 
els relevant for a unified description of diffusion, second sound, and ballistic 
regimes of heat propagation are intensively discussed (see [27, 28] and ref- 
erences therein). To be specific, recall the simplest and typical model of 
the phonon transport [27]. Let e(x,t) and p(x, t) be small deviations of 
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the energy density and of the energy flux of the phonon field from their 
equilibrium values, respectively. Then 

d t e = -c 2 V-p, (125) 
dtp = -lve-—p. (126) 

O Tr 

Here c is the Debye velocity of phonons, and tr is the characteristic time 
of resistive processes. Equations (125) can be derived from the Boltzmann- 
Peierls kinetic equation, within the relaxation time approximation, by a 
method similar to the Grad method [27]. Eqs. (125) provide the simplest 
model of a coupling between the hydrodynamic variable e and the non- 
hydrodynamic variable p, allowing for a qualitative description of both the 
diffusion and the second sound. Following the standard argumentation [27], 
we observe the two limiting cases: 1). As tr — > 0, the equation (126) yields 
the Fourier relation p = — ^trVc which closes the equation (125) to give 
the diffusion equation: 

dte + ^T R c 2 Ae = 0. (127) 

2). As tr — ► oo, the equation (126) yields dtp = — jVe, and the equation 
(125) closes to give the wave equation: 

d 2 t e + \c 2 Ae = 0. (128) 

The equation (127) describes the usual diffusive regime of the heat propaga- 
tion, while the equation (128) is relevant to the (undamped) second sound 
regime with the velocity U2 = c/y/3, and are both closed with respect to the 
variable e. 

However, even within the simplest model (125), the problem of closure 
remains unsolved in a systematic way when tr is finite. The natural way 
of doing so is provided by the Chapman-Enskog method. In the situation 
under consideration, the Chapman-Enskog method yields an extension of 
the diffusive transport to finite values of the parameter tr, and leads to an 
expansion of the non-hydrodynamic variable p in terms of the hydrodynamic 
variable e. With this, if we are able to make this extension of the diffusive 
mode exactly, we could learn more about the transition between the diffusion 
and second sound (within the frames of the model). 
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The Chapman-Enskog method, as applied to the model (125), results in 
the following series representation: 

oo 

p=Ep w . ( 129 ) 

n=0 

where the coefficients p^ are due to the Chapman-Enskog recurrence pro- 
cedure, 

n-l 

p^ = -r R Y,d { r ] P {n -^ m \ (130) 

m=0 

while the Chapman-Enskog operators d\ m ^ act on e as follows: 

5 M e = _ c 2 v . p M (m) 

Finally, the zero order term reads: p^ = — |trVc, and leads to the Fourier 
approximation of the energy flux. 

To sum up the series (129) in a closed form, we, will specify the non- 
linearity appearing in equations (130) and (131). The coefficients 

p(n) in 

equations (129) and (130) have the following explicit structure for arbitrary 
order n > 0: 

p {n) = a n A n Ve, (132) 

where the real- valued and yet unknown coefficients due to the recur- 

rence procedure (130), and (131). Indeed, the form (132) is true for n = 
(do = — gTfl). Let us assume that (132) is proven up to the order n — 1. 
Then, computing the nth order coefficient p^ , we derive: 

n-l 

P (n) = -TR^dt^n-l-rrA^-^Ve (133) 
m=0 

n-l 

= -tr ]T a„_ 1 _ m A (ri - 1 " m) V (-c 2 a m V • VA m e) 

m=0 

= t r c 2 I a n-i- m a m > A n Ve. 

lm=0 J 

The last expression has the form (132). Thus, the Chapman-Enskog pro- 
cedure for the model (125) is equivalent to the following nonlinear recurrence 
relation in terms of the coefficients a n : 
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n-1 

a n = TRC 2 £ "n-l-mflm, (134) 
m=0 

subject to the initial condition ao = — \t~r- Further, it is convenient to make 
the Fourier transform. Using p = p k ex.p{ik-x} and e = e k exp{ik-x}, where 

k is the real-valued wave vector, we derive in (132): pj^ = a n ik(—k 2 ) n e k , 
and 

p k = ikA(k 2 )e k , (135) 

where 

oo 

A(k 2 ) = Y / a n (-k 2 ) n . (136) 

n=0 

Thus, the the Chapman-Enskog solution (129) amounts to finding the func- 
tion A{k 2 ) represented by the power series (136). If the function A is known, 
the exact Chapman-Enskog closure of the system (125) amounts to the fol- 
lowing dispersion relation of plane waves ~ exp{oj k t + ik ■ x}: 

u k = c 2 k 2 A(k 2 ). (137) 

Here uj k is a complex- valued function of the real- valued vector k: Reu) k is 
the attenuation rate, Imuj k is the frequency. 

Multiplying both the equations in (134) with (—k 2 ) n , and performing a 
summation in n from 1 to infinity, we get: 

oo n 

A-a = -r R c 2 k 2 ]T ]T a n - m {-k 2 ) n - m a m {-k 2 ) m , 

n=0 m=0 

Now we notice that 

TV n 

Jim £ £ a n . m (-k 2 ) n - m a m (-k 2 r = A 2 , 

n=U m=U 

Accounting ao = — ^tr, we come to a quadratic equation for the function 
A: 

T R c 2 k 2 A 2 + A + ^t r = 0. (138) 

Further, a selection procedure is required to choose the relevant root of the 
equation (138). Firstly, recall that all the coefficients a n (132) are real- 
valued by the sense of the Chapman-Enskog method (130) and (131), hence 
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the function A (136) is real-valued. Therefore, only the real-valued root of 
the equation (138) is relevant to the Chapman-Enskog solution. The first 
observation is that the equation (138) has no real- valued solutions as soon 
as k is bigger than the critical value k c , where 



Secondly, there are two real-valued solutions to the equation (138) at k < 
k c . However, only one of them satisfies the Chapman-Enskog asymptotic 
lim k ^ A(k 2 ) = -\tr. 

With the two remarks just given, we finally derive the following exact 
Chapman-Enskog dispersion relation (137): 



The Chapman-Enskog dispersion relation corresponds to the extended 
diffusion transport, and it comes back to the standard Fourier approximation 
in the limit of long waves k/k c <C 1. The Chapman-Enskog solution does 
not exist as soon as k/k c > 1. In the point k c , the extended diffusion branch 
crosses one of the non-hydrodynamic branches of the (125). For larger k, the 
extended diffusion mode and the critical non-hydrodynamic mode produce 
a pair of complex conjugated solutions with the real part equal to — 
The imaginary part of this extension after k c has the asymptotics ±iu 2 k, as 
k — > oo, and where 112 = c/s/3 is the (undamped) second sound velocity in 
the model (125) (see equation (128)). Though the spectrum of the original 
eq. (125) continues indeed after k c , the Chapman-Enskog method does not 
recognize this extension as a part of the hydrodynamic branch, while the 
second sound regime is born from the extended diffusion after the coupling 
to the critical non-hydrodynamic mode. 

Finally, let us consider the opportunities provided by the Newton method 
as applied to the invariance equation. First, the invariance equation can be 
easily obtained in a closed form here. Consider again the expression for 
the heat flux in terms of the energy density (135), p k = ikA(k 2 )ek, where 
now the function A is not thought as the Chapman-Enskog series (136). 
The invariance equation is a constraint on the function A, expressing the 
form- invariance of the heat flux (135) under both the dynamic equations 
(125) and (126). Computing the time derivative of the function (135) due 





(140) 



none k > k c 
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to equation (125), we obtain: 



5 macro pfc = ikA ( k ^d t e k = c 2 k 2 A 2 ike k . (141) 

On the other hand, computing the time derivative of the same function due 
to equation (126), we have: 

df CT °Pk = ~\ike k - -Aike k . (142) 

Equating the expressions (141) and (142), we come to the desired invari- 
ance equation for the function A. This equation coincides with the exact 
Chapman-Enskog equation (138). 

As the second step, let us apply the Newton method to the invariance 
equation (138), taking the Euler approximation (^4^ = 0) for the initial 
condition. Rewriting the equation (138) in the form F(A, k 2 ) = 0, we come 
to the following Newton iterations: 

m 1 h " ' (A n+ i - An) + F(A n , k 2 ) = 0. (143) 



dA 



A — A n 



The first two iterations give: 



TrA 1 = -i (144) 

The first Newton iteration (144) coincides with the first term of the 
Chapman-Enskog expansion. The second Newton iteration (145) is a ratio- 
nal function with the Taylor expansion coinciding with the Chapman-Enskog 
solution up to the super-Burnett term, and it has a pole at yi = V2. The 
further Newton iterations are also rational functions with the relevant poles 
in the points y n , and the sequence of this points tends very rapidly to the 
location of the actual singularity y c = 1 (^3 w 1.17, 1/4 w 1.01, etc.). 



3.6.2 Inclusion of normal processes 

The account for normal processes in frames of the semi-hydrodynamical 
models [27] leads to the following generalization of the Eq. 125 (written in 
Fourier variables, in the one-dimensional case): 
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dtek = -ikc 2 p k , (146) 

dtPk = ~-ike k -ikN k p k , (147) 

3 t r 

8 t N k = -^-ikc 2 Pk --N k . (148) 



Here r = tntr/(tn + tr), tjv is the characteristic time of normal pro- 
cesses, and N k is the additional field variable. Following the principle of 
invariance as explained in the preceding section, we write the closure rela- 
tion for the non-hydrodynamic variables p k and N k as: 

p k = ikA k e k , N k = B k e k , (149) 

where A k and Bk are two unknown functions of the wave vector k. Further, 
following the principle of invariance as explained in the preceding section, 
each of the relations (149) should be invariant under the dynamics due to 
Eq. (146), and due to Eqs. (147) and (148). This results in two equations 
for the functions Ak and B k : 

k 2 c 2 A 2 k = -L Ak -B h -\, (150) 

TR 6 

k 2 c 2 A k B k = --B k + ^-k 2 c 2 A k . 

T 15 

When the energy balance equation (146) is closed with the relation (149), 
this amounts to a dispersion relation for the extended diffusion mode, to k = 
k 2 (?A k , where A k is the solution to the invariance equations (150), subject 
to the condition A k — > as k — > 0. Resolving equations (150) with respect 
to A k , and introducing A k = k 2 c 2 Ak, we arrive at the following: 

= + _ _ i ^ (151) 

5 + 9rAk o 

The invariance equation (151) is completely analogous to the Eq. (138). 
Written in the form (151), it allows for a direct investigation of the critical 
points. For this purpose, we find zeroes of the derivative, d^(A k )/dA k = 0. 
When the roots of the latter equation, A c k , are found, the critical values of the 
wave vector are given as — (l/3)/c^c 2 = <&(A° k ). The condition d^{A k )/dAk = 
reads: 

18r 2 r^ + 3r(3r + 8t r )A^ + 10(r + r R )A k + 5 = 0. (152) 
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Let us consider the particularly interesting case, e = tn/tr <C 1 (the 
normal events are less frequent than resistive). Then the real- valued root 
of the equation (152), A k (e), corresponds to the coupling of the extended 
diffusion mode to the critical non-hydrodynamic mode. The corresponding 
modification of the critical wave vector k c (139) due to the normal processes 
amounts to a shifts towards shorter waves, and we derive: 

lkc(e)] 2 = k 2 + (153) 
3.6.3 Account for anisotropy 

The above examples concerned the isotropic Debye model. Let us consider 
the simplest anisotropic model of a cubic media with a longitudinal (L) and 
two degenerated transverse (T) phonon modes, taking into account resistive 
processes only. Introduce the Fourier variables, e k , el, pi, and p%, where 
e k = e£ +2ejT is the Fourier transform of the total energy of the three phonon 
modes (the only conserving quantity), while the rest of variables are specific 
quantities. The isotropic model (125) generalizes to give [27]: 

d t e k = -ic 2 L k ■ pi - 2%c\k ■ pi, (154) 
3 t el = -ic 2 T k -p T k+\ [c£(e* - 2e T k ) ~ 44] , (155) 

3 tV L k = -Uk(e k - 2el) - \pt (156) 

<3 Tr- 



<R 



dtPl = -giMfc - -^rPk , (157) 



where A = t^Cj. + 1t^c\. The term containing the factor A -1 corresponds 
to the energy exchange between the L and T phonon modes. The invariance 
constraint for the closure relations, 

Pk = ikA k e k , pi = ikB k e k , el = X k e k , (158) 

result in the following invariance equations for the fc-dependent functions 
A k , B k , and X k : 



k 2 c 2 L A 2 k + 2k 2 c 2 T A k B k = -l-A k -\(l-2X k ), (159) 
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2k 2 c 2 T B 2 + k 2 c L B k A k = -—B k --X k , (160) 

T R 6 

X k (k 2 c 2 L A k + 2k 2 clB k ) = c\k 2 B k + - [c\ - X k (2c| + 4)] .(161) 

When the energy balance equation (154) is closed with the relations (158), 
this leads to the dispersion relation for the extended diffusion mode, uj k = 
A k + 2B k , where the functions A k = k 2 c 2 L A k , and B k = k 2 c^B k , satisfy 
the condition: A k — > 0, and B k — ► 0, as k — > 0. The resulting dispersion 
relation is rather complicated in the general case of the four parameters of 
the problem, cl, ct, t r and t r . Therefore, introducing a function Y k = 
A k + 2B k , let us consider the following specific situations of closed equations 
for the Y k on the basis of the invariance equations (159): 

(i). cl = ct = c, t r = t r = t r (complete degeneration of the parame- 
ters of the L and T subsystems): The system (159) results in two decoupled 
equations: 

Y k (r R Y k + l) + \k 2 c 2 T R = 0, (162) 



(r R Y k 



3 

+ l) 2 + h 2 c 2 r 2 R = 0. (163) 



The equation (162) coincides with the eq. (138) for the isotropic case, 
and its solution defines the coupling of the extended diffusion to a non- 
hydrodynamic mode. The equation (163) does not have a solution with the 
required asymptotic Y k — > as k — > 0, and is therefore irrelevant to the fea- 
tures of the diffusion mode in this completely degenerated case. It describes 
the two further propagating and damped non-hydrodynamic modes of the 
Eqs. (154). The nature of these modes, as well of the mode which couples 
to the diffusion mode well be seen below. 

(ii). cl = ct = c, ^ t r (nondegenerate characteristic time of resistive 
processes in the L and the T subsystems): 



2 ,2 



Y k (4Y k + l) + \k 2 c 2 r^ x \^r' R Y k + 3) (r T R Y k + l) + h 2 c 2 r^r' R 

J (164) 

where t' r = 2t r +t r . As t r —t r — > 0, the Eq. (164) tends to the degenerated 
case (162). At k = 0, t r ^ t r , there are four solutions to the Eq. (164). 
The yo = is the hydrodynamic solution indicating the beginning of the 
diffusion mode. The two non-hydrodynamic solutions, Yq = —1/t r , and 
Yq = —1/t r , Yq = —3/t' r , are associated with the longitudinal and the 



52 



transverse phonons, respectively. The difference in relaxational times makes 
the latter transverse root nondegenerate, instead there appears the third 
non-hydrodynamic mode, Yq = —3/t' r . 

(iii). cl / ct, tJi = Tft = tr (nondegenerate speed of the L and the T 
sound). 

Y k (T R Y k + l) + l -k 2 c 2 L T R 

"(165) 

As ct — cl —>■ 0, the Eq. (165) tends to the degenerated case (162). How- 
ever, this time the non-hydrodynamic mode associated with the transverse 
phonons is degenerated at k = 0. 

Thus, we are able to identify the modes in the equations (162) and (163). 
The non-hydrodynamic mode which couples to the extended diffusion mode 
is associated with the longitudinal phonons, and is the case (162). The case 
(163) is due to the transverse phonons. In the nondegenerate cases, (164) 
and (165), the both pairs of modes become propagating after certain critical 
values of k, and the behavior of the extended diffusion mode is influenced 
by all the three non-hydrodynamic modes just mentioned. It should be 
stressed, however, that second sound mode which is the continuation of the 
diffusion mode [24], is due to the eq. (162). The results of the above analysis 
lead to the following discussion: 

(i) The examples considered above indicate an interesting mechanism of 
a kinetic formation of the second sound regime from the extended diffusion 
with the participation of the non-hydrodynamic mode. The onset of the 
propagating mode shows up as the critical point of the extension of the 
hydrodynamic solution into the domain of finite k, which was found within 
the Chapman-Enskog and equivalent approaches. This results concern the 
situation at the high-temperature edge of the Guyer-Krumhansl window, 
and are complementary to the coupling between the transversal ballistic 
mode and second sound at the low-temperature edge [29]. 

(ii) The crossover from the diffusion-like to the wave-like propagation 
was previously found in [30] in frames of exact Chapman-Enskog solution 
to the Boltzmann equation for the Lorentz gas model [7], and for simi- 
lar models of phonon scattering in anisotropic disordered media [31]. The 
characteristic common feature of the models studied in [7, 30, 31] and the 
models [27] is the existence of a gap between the hydrodynamic (diffusive) 
and the non-hydrodynamic components of the spectrum. Therefore, one 
can expect that the destruction of the extended diffusion is solely due to 



(r R Y k + l) 



2 + ^ 2 c 2 T r 2 R 



2 ,< 



TR 



'L\ T 



2c? + cl 



(r R Y k + l) = 0. 
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the existence of this gap. In applications to the phonon kinetic theory this 
amounts to the introduction of the relaxation time approximation. In other 
words, we may expect that the mechanism of crossover from diffusion to sec- 
ond sound in the simple models [27] is identical to what could be found from 
the phonon-Boltzmann kinetic equation within the relaxation time approx- 
imation. However, a remark is in order since the original (i. e., without the 
relaxation time approximation) phonon kinetic equations are gapless (cf., 
e. g., [26]). On the other hand, most of the works on heat propagation in 
solids do explore the idea of the gap, since it is only possible to speak of 
the diffusion if such a gap exists. To conclude this point, the following gen- 
eral hypothesis can be expressed: the existence of the diffusion ( and hence 
of the gap in the relaxational spectrum) leads to its destruction through the 
coupling to a non-hydrodynamic mode. 

3.7 Nonlinear Grad equations 

In the preceding sections, the Chapman-Enskog and other methods were 
probed explicitly for the linearized Grad equations far beyond the usual 
Navier-Stokes approximation. This was possible, first of all, because the 
problem of the reduced description was shaped into a rather simple algebraic 
form. Indeed, the algebraic structure of the stress tensor ak{pk,Uk,Tk,k) 
and of the heat flux q^iPk, Uf~, T^, k) was fairly simple. However, when we 
attempt to extend the approach onto the nonlinear Grad equations, the 
algebraic structure of the problem is no longer simple. Indeed, when we 
proceed along the lines of the Chapman-Enskog method, for example, the 
number of types of terms, Vm, VVu, (Vm) 2 , VTVp, and so on, in the 
Chapman-Enskog coefficients cr^ and demonstrates the combinatorial 
growth with the order n. 

Still, a progress is possible if we impose some rules of selection of rele- 
vant terms. As applied to the Chapman-Enskog expansion, these selection 
rules prescribe to keep only the contributions coming from the terms with a 
definite structure in each order and q^ n \ and to ignore all other terms. 
This approach can be linked again with the partial summation rules for the 
perturbation series in many-body theories, where usually terms with a def- 
inite structure are summed instead of a series in a whole. Our viewpoint 
on the problem of extension of the hydrodynamics in the nonlinear case 
can be expressed as follows: The exact extension seems to be impossible, 
and, moreover, quite useless because of the lack of a physical transparency. 
Instead, certain sub-series of the Chapman-Enskog expansion, selected on 
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clear physical grounds, may lead to less complicated equations, which, at 
the same time, provide an extension of a certain subclass of hydrodynamic 
phenomena. This viewpoint is illustrated in this section by considering a 
sub-series of the Chapman-Enskog expansion which gives the dominating 
contribution when the flow velocity becomes very large (and thus it is rel- 
evant to a high-speed subclass of hydrodynamic phenomena such as strong 
shock waves). 

The approach to the Chapman-Enskog series for the nonlinear Grad 
equations just mentioned, and which was based on a diagrammatic repre- 
sentation of the Chapman-Enskog method, has been attempted earlier in 
[14]. In this section, however, we will take the route of the dynamic invari- 
ance equations which leads to the same results more directly. 

3.7.1 The dynamic viscosity factor 

The starting point is the set of one-dimensional nonlinear Grad equations for 
the hydrodynamic variables p, u and T, coupled to the non-hydrodynamic 
variable a, where a is the xx-component of the stress tensor: 

dtp = -d x (pu); (166) 

d t u = -ud x u - p~ x d x p - p~ x d x o\ (167) 

d t T = -ud x T - (2/3)Td x u- (2/3)p~ 1 (jd x u; (168) 

d t a = -ud x a-(4/3)pd x u-(7/3)ad x u-jj^a. (169) 

Here p(T) is the temperature-dependent viscosity coefficient. We will adopt 
the form p(T) = aT 7 , which is characteristic to the point-center models of 
particle's collisions, where 7 varies from 7 = 1 (Maxwellian molecules) to 
7 = 1/2 (hard spheres), and where a is a dimensional factor. 

Even in this model, the Chapman-Enskog expansion appears to be ex- 
ceedingly complicated in the full setting. Therefore, we address another, 
more restrictive problem: What is the leading correction to the Navier- 
Stokes approximation when the characteristic value of the average velocity 
is comparable to the thermal velocity? 

Our goal is to compute the correction to the Navier-Stokes closure re- 
lation, o"ns = ~~ (4/3)pd x u, for high values of the average velocity. Let us 
consider first the Burnett correction as from the Eq. (166): 

<tb = -\pd x u + 8(2 ~ 7 V p- 1 ($r") 2 " \v 2 p- l d x {p~ l d xP ). (170) 
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The correction of the desired type is given by the nonlinear term propor- 
tional to (d x u) 2 . Each further nth term of the Chapman-Enskog expansion 
contributes, among others, the nonlinear term proportional to (d x u) n+1 . 
Such terms can be named high-speed terms since they dominate the rest 
of the contributions in each order of the Chapman-Enskog expansion when 
the characteristic average velocity is comparable to the heat velocity. In- 
deed, if U is a characteristic mean velocity, and u = Uu, where u is di- 
mensionless, then the term (d x u) n+1 receives the factor JJ n+l which is the 
highest possible order of U among the terms available in the nth order of 
the Chapman-Enskog expansion. Simple dimensional analysis leads to a 
conclusion that such terms are of the form [ig n d x u, where g = p~ 1 [id x u is 
dimensionless. Therefore, the Chapman-Enskog expansion for the function 
a may be formally rewritten as: 



At 



-^-?lg + r 2 g 2 + . . . + r n g n + . . )d x u + . . . (171) 



The series in the brackets is the collection of the high-speed contributions 
of interest, coming from all the orders of the Chapman-Enskog expansion, 
while the dots outside the brackets stand for the terms of other nature. Thus, 
the high-speed correction to the Navier-Stokes closure relation in frames of 
the Grad equations (166) takes the form: 

<j n i = -vR{g)d x u, (172) 

where R(g) is the function represented by a formal subsequence of Chapman- 
Enskog terms in the expansion (171). The function R can be viewed also as 
a dynamic modification of the viscosity u due to the gradient of the average 
velocity. 

We will now turn to the problem of an explicit derivation of the function 
R (172). Following the principle of dynamic invariance, we, first, compute 
the microscopic derivative of the function a n \ by substituting the expression 
(172) into the right hand side of the Eq. (169): 

4 7 p (47 1 

dl mcro a ld = -ud x a n i--pd x u--a n id x u-—^a n i = j-- + -gR + R j pd x u+. . . , 

(173) 

where dots denote the terms irrelevant to the closure relation (172) (such 
terms appear, of course, because (172) is not the exact closure relation). 
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Second, computing the macroscopic derivative of the closure relation 
(172) due to the equations (166), (167), and (168), we have: 

JO 

a t macr V nl = -[d t n(T)]Rd x u - n(T) — [d t g]d x u - fi(T)Rd x [d t u\. (174) 

In the latter expression, the time derivatives of the hydrodynamic variables 
should be replaced with the right hand sides of the Eqs. (166), (167), and 
(168), where, in turn, the function a should be replaced by the function a n \ 
(172). After some computation, we come to the following: 



d t macr V nI = {gR+ |(1 - gR) X (jgR+^-l)g2^j}pd xU +... (175) 

Here dots stand again for the terms irrelevant to the present analysis. 

Equating the relevant terms in Eqs. (173) and (175), we come to the 
invariance equation for the function R: 



(l- 7 ) 5 2 (1-^)^+7^ + 



R-2 = 0. (176) 



For Maxwell molecules (7 = 1), Eq. (176) simplifies considerably, and is 
algebraic: 

g 2 R 2 + (l+g]R- 2 = 0. (177) 



2 

The solution which recovers the Navier-Stokes closure relation in the limit 
of small g reads: 

flMM = zi^±M^W gW. (178) 

Function Rmm (178) is demonstrated in Fig. 10. Notice that the function 
Rmm is positive for all values of its argument g, as is appropriate for the vis- 
cosity factor, while the Burnett approximation to the function Rmm violates 
the positivity. 

For other models (7 / 1), the invariance equation (176) is a rather com- 
plicated nonlinear ODE with the initial condition R(0) = 4/3 (the Navier- 
Stokes condition). Several ways to derive analytic results are possible. One 
possibility is to expand the function R into powers of g, in the point g = 0. 
This way will bring us back to the original sub-series of the Chapman-Enskog 
expansion (see Eq. (171)). Instead, we take the opportunity offered by the 
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parameter 7. Introduce another parameter (3 = 1 — 7, and consider the 
expansion: 

R(P,g) = R (g) + PRi(g) + 2 R 2 (g) + ■■■■ 

Substituting this expansion into the invariance equation (176), we derive 
Ro(g) = Rmm {g) j 

o ^ (l-gRo)gRo + (3/2)g 2 (dR /dg) 

Rl{9) = w^Tg-rm — ' (179) 

etc. That is, the solution for other models is constructed in a form of 
a series with the exact solution for the Maxwell molecules as the leading 
term. For hard spheres (/? = 1/2), the result to the first-order term reads: 
i? HS « Rmm + (l/2)i?i- The resulting approximate viscosity factor is given 
in Fig. 11. The features of the approximation obtained are qualitatively the 
same as in the case of the Maxwell molecules. 



3.7.2 Attraction to the invariant set 

Above, we have derived a correction to the Navier-Stokes expression for the 
stress <t, in the one-dimensional case, for large values of the average velocity 
u. This correction has the form a = —fiR(g)d x u, where g oc d x u is the 
longitudinal rate. The viscosity factor R(g) is a solution to the differential 
equation (176), subject to a certain initial condition. Uribe and Piha [32] 
have indicated some interesting features of the invariance equation (176) for 
the model of hard spheres. In particular, they have found that a numerical 
integration from the initial point into the domain of negative longitudinal 
rates is very difficult. What happens to the relevant solution at negative 
values of gl 

Let us denote as P = (g, R) the points of the (g, R) plane. The relevant 
solution R(g) emerges from the point Pq = (0, 4/3), and can be uniquely con- 
tinued to arbitrary values of g, positive and negative. This solution can be 
constructed, for example, with the Taylor expansion, and which is identical 
with the relevant sub-series of the Chapman-Enskog expansion. However, 
the difficulty in constructing this solution numerically for g < originates 
from the fact that the same point Pq is the point of essential singularity 
of other (irrelevant) solutions to the Eq. (176). Indeed, for \g\ <C 1, let us 
consider R(g) = R(g) + A, where R(g) = (4/3) + (8/9)(j-2)g is the relevant 
solution for small \g\, and A (5) is a deviation. Neglecting in the Eq. (176) 
all regular terms (of the order g 2 ), and also neglecting g A in comparison to 
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A, we derive the following equation: (1 — r ))g 2 {dA/dg) = — (3/2)A. The 
solution is A(g) = A(<?o) exp[a((/ _1 — g^ 1 )], where a = (3/2)(l — 7) -1 . The 
essential singularity at g = is apparent from this solution, unless A(<?o) / 
(that is, no singularity exists only for the relevant solution, R = R). Let 
A(<7q) 7^ 0. If g < 0, then A — > 0, together with all its derivatives, as <? — ► 0. 
If g > 0, the solution blows up, as # — ► 0. 

The complete picture for 7 ^ 1 is as follows: The lines g = and P = 
(<7,<7 _1 ) define the boundaries of the basin of attraction A = ^-IJ^+j where 
A- = {P\ - 00 < g < 0,R > g- 1 }, and A + = {P\oo > g > 0, R < g' 1 }. 
The graph G = (g, R(g)) of the relevant solution belongs to the closure of A, 
and goes through the points Pq = (0, 4/3), P_ = (—00, 0), and P + = (00, 0). 
These points at the boundaries of A are the points of essential singularity of 
any other (irrelevant) solution with the initial condition P G A, P £ Af] G. 
Namely, if P G A + , P ^ A + f] G, the solution blows up at Pq, and attracts 
to P + . If P G A-, P A^ H G, the solution blows up at P_, and attracts 
to P . 

The above consideration is supported by a numerical study of the Eq. 
(176). In Fig. 12, it is demonstrated how the dynamic viscosity factor R(g) 
emerges as the attractor of various solutions to the invariance equation (176) 
[the case considered corresponds to hard spheres, 7 = 1/2]. The analytical 
approximation (179) is also shown in Fig. 12 for the sake of comparison. It 
provides a reasonable global approximation to the attractor for both positive 
and negative g. We conclude with a discussion. 

(i) The main feature of the above example of extending the hydrody- 
namic description into a highly non-equilibrium and nonlinear domain can 
be expressed as follows: this is an exact partial summation of the Chapman- 
Enskog expansion. "Partial" means that the relevant high-speed terms, 
dominating the other contributions in the limit of the high average velocity, 
were accounted to all the orders of the original Chapman-Enskog expansion. 
"Exact" means that, though we have used the formally different route, the 
result is indeed the sum of the relevant sub-series of the original Chapman- 
Enskog expansion. In other words, if we now expand the function -Rmm(<?) 
(178) in powers of g, in the point g = 0, we come back to the corresponding 
series inside the brackets in the Eq. (171). That this is indeed true can be 
checked up to the few lower orders straightforwardly, though the complete 
proof requires a more involved analysis, and will be reported elsewhere. As 
the final comment to this point, we would like to stress a certain similarity 
between the problem considered above and the frequent situations in many- 
body problems: there is no a single leading term but instead there is the 
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leading sub-series of the perturbation expansions, under certain conditions. 

(ii) Let us discuss briefly the features of the resulting hydrodynamics. 
The hydrodynamic equations are now given by Eqs. (166), (167), and (168), 
where a is replaced with a n \ (172). First, the correction concerns the non- 
linear regime, and thus the linearized form of new equations coincides with 
the linearized Navier-Stokes equations. Second, the solution (178) for the 
Maxwellian molecules and the result of approximation (179) for other mod- 
els suggests that the modified viscosity gives a vanishing contribution 
in the limit of very high values of the average velocity. This feature seems 
to be of no surprise: if the average velocity is very high in comparison to 
other characteristic velocities (in our case, to the heat velocity), no mech- 
anisms of momentum transfer are relevant except for the transfer with the 
stream. However, a cautious remark is in order since the original "kinetic" 
description are Grad equations (166) and not the Boltzmann equation. 

(iii) The invariance equation (176) defines the relevant physical solution 
to the viscosity factor for all values of g, and demonstrates an interesting 
phase-space behavior similar to those of finite-dimensional dynamical sys- 
tems. 

4 Conclusions 

Up to now, the problem of the exact relationship between kinetics and hy- 
drodynamics remains unsolved. All the methods used to establish this re- 
lationship are not rigorous, and involve approximations. In this work, we 
have considered situations where hydrodynamics is the exact consequence 
of kinetics, and in that respect, a new class of exactly solvable models of 
statistical physics has been established. 

The main lesson we can learn from the exact solution is the following: 
The Chapman-Enskog method is the Taylor series expansion approach to 
solving the appropriate invariance equation. Alternative iteration methods 
are much more robust for solving this equation. Therefore, it seems quite 
important to develop approaches to the problem of reduced description based 
on the principle of dynamic invariance rather than on particular methods of 
solving the invariance equations. The exact solutions where these questions 
can be answered in all the quantitative details provide a sound motivation 
for such developments. Interested reader is directed to the paper [23] where 
a general formulation of the reduced description from the standpoint of 
constructing invariant manifolds has been developed, as well as to some 
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recent applications in the problems of Boltzmann kinetic theory, chemical 
kinetics, and polymer kinetic theory are given [33, 34, 35, 36]. 
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Figure 1: Real- valued root of Eq. (33) as a function of k 2 
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Figure 2: Attenuation rates for the IDIOM Grad system. Solid: Exact 
summation of the Chapman-Enskog Expansion. Dots: The Navier-Stokes 
approximation. Dash: The super-Burnett approximation. Circles: Hydro- 
dynamic and non-hydrodynamic modes of the IDIOM Grad system. 
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Figure 3: Attenuation rates for the 3-D10M Grad system. Bold: The acous- 
tic branch, exact summation. Dots: The acoustic branch, Navier-Stokes 
approximation. Circles: The acoustic branch, super-Burnett approxima- 
tion. Solid: The diffusion branch, exact summation. Dash: The critical 
mode of the 3D10M Grad system. 
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Figure 4: Attenuation rates for the partial summing. Solid: The regularized 
Burnett approximation. Dash: The regularized super-Burnett approxima- 
tion. Circles: The super-Burnett approximation. Dots: The exact summa- 
tion. 
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Figure 5: Attenuation rates for the Newton method with the Navier-Stokes 
approximation as the initial condition. Dots: The Navier-Stokes approxi- 
mation. Solid: The first and the second iterations of the invariance equa- 
tion. Circles: The exact solution to the invariance equation. Diamond: The 
super-Burnett approximation. 
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Figure 6: Attenuation rates with the regularized Burnett approximation as 
the initial condition for the Newton method. Dots: The regularized Burnett 
approximation, or the first Newton iteration with the Euler initial condition 
(see text). Solid: The first and the second Newton iterations with the 
regularized Burnett approximation as the initial condition. Circles: The 
exact solution to the invariance equation. 
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Figure 7: Attenuation rates with the regularized super-Burnett approxima- 
tion as the initial condition for the Newton method. Dots: The regularized 
super-Burnett approximation. Solid: The first and the second Newton iter- 
ations. Circles: The exact solution to the invariance equation. 
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Figure 8: The diffusion mode with the Euler initial approximation for the 
invariance equation. Dot-dash: The the first iteration. Dots: The second 
iteration. Circles: The third iteration. Solid: The exact solution. Dash: 
The critical mode. 
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Figure 9: Solutions to the Eq. (123). Circles: Numerical solutions. Dots: 
The first Newton iteration. Solid: The 4th Newton iteration. 
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Figure 10: Viscosity factor R(g) for Maxwell molecules. Solid: exact solu- 
tion. Dash: the Burnett approximation. Dots: the Navier-Stokes approxi- 
mation. 
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Figure 11: Viscosity factor R{g) for hard spheres. Solid: the first approxi- 
mation. Dash: exact solution for Maxwell molecules. 
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Figure 12: Viscosity factor as an attractor. Solid lines: numerical integration 
with various initial points (crosses). Two poorly resolved lines correspond 
to the initial conditions (—100,0) and (—100,3). Circles: Taylor expan- 
sion to the 5th order. Dots: the analytical approximation of (179). Dash: 
boundaries of the basin of attraction. 
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